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In a number of applications, particularly in financial and actuar- 
ial mathematics, it is of interest to characterize the tail distribution of 
a random variable V satisfying the distributional equation V — f{V), 



'~~i for some random function /. This paper is concerned with computa- 

r — tional methods for evaluating these tail probabilities. We introduce 

T-H a novel dynamic importance sampling algorithm, involving an expo- 
nential shift over a random time interval, for estimating such rare 

' J ' event probabilities. We prove that the proposed estimator is: (i) con- 

Mh sistent; (ii) strongly efficient; and (iii) optimal within a wide class of 

Mh dynamic importance sampling estimators. Moreover, using extensions 

ri of ideas from nonlinear renewal theory, we provide a precise descrip- 

• *— * tion of the running time of our algorithm. To establish these results, 

_i we develop new techniques concerning the convergence of moments 

S of stopped perpetuity sequences, and the first entrance and last exit 

times of associated Markov chains on R. We illustrate our methods 

I with a variety of numerical examples that demonstrate the ease and 

►^ scope of the implementation. 

00 

C^ 1. Introduction. Recently there has been much attention focused on the 

study of stochastic fixed point equations. This interest has been largely fu- 
eled by contemporary applications, including aspects of financial time series 
modeling, actuarial mathematics, the study of perpetuities, and the anal- 

^~~i ysis of algorithms. An original motivation for this article came from one 

J> such application, namely the desire to provide a useful and implementable 

*k> algorithm — that is theoretically justified — for calculating the probability of 
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ruin for an insurance company under the Cramer-Lundberg model, but in 
the presence of stochastic investments. A cfosely related problem is that of 
obtaining precise estimates for the tail probabilities in the GARCH(1,1) and 
ARCH(l) financial time series models, which are needed for Value-at-Risk or 
expected shortfall calculations in risk management. In spite of the relevance 
of these and a number of related problems, effective computational methods 
for evaluating these rare event probabilities have, to date, been investigated 
rather poorly. 

This paper is concerned with importance sampling methods for calculating 
the stationary tail probabilities of a stochastic fixed point equation (SFPE). 
In general, an SFPE assumes the form 

(1.1) V^fiV), 

where / is a random function satisfying certain regularity conditions and is 
independent of V. When f{v) = Av-\-B, where E [log A] < 0, this problem has 
a long history beginning with the works of Solomon (1972), Kesten (1973), 
and Vervaat (1979). In this work, we consider a generalization of their model, 
namely Letac's (1986) Model E, whose SFPE is given by 

(1.2) V = Amax{D,V)+B. 

This recurrence equation is known to be quite general in the sense that it 
characterizes a wide variety of mathematical problems of applied impor- 
tance. These applications include the ruin problem with investments, the 
GARCH(1,1) and ARCH(l) financial time series models, the AR(1) time se- 
ries model with random coefficients, and the classical GI/G/1 queue, among 
others. 

In a recent work, Collamore and Vidyashankar (2011) extended earlier 
results of Kesten (1973) and Goldie (1991), establishing 

(1.3) lim u^P {V > u) = C 

for finite positive constants C and S,, where C is identified as the difference of 
the ^th moments of a perpetuity sequence and another "backward sequence," 
and ^ is identified as the nonzero solution to the equation E [A^] = 1. To 
obtain an exact estimate for P {V > u), however, it is natural to turn to 
computational methods and particularly to importance sampling. 

In large deviation problems for sums of i.i.d. or Markov random walks, 
importance sampling methods have been developed by a number of authors. 
Following the lines of Hammersley and Handscomb (1964), an early work 
is Siegmund (1976), who introduced an algorithm for computing the level 



IMPORTANCE SAMPLING FOR SFPE 3 

crossing probabilities in the sequential probability ratio test. This method was 
later generalized to other large deviation problems in, e.g., Asmussen (1985), 
Lehtonen and Nyrhinen (1992), Bucklew et al. (1990), CoUamore (2002), Chan 
and Lai (2007), and Blanchet and Liu (2010). 

The problem we consider here is different, since we study processes with 
both multiplicative and additive components. Indeed, if one directly adopts 
importance sampling methods for either light- or heavy-tailed sums, then 
the resulting estimator will not be efficient. Instead, we propose a different 
approach, involving a dual change of measure over a random time interval, 
simulated over the finite-time excursions of {Vn} emanating from a certain set 
C C M and returning to this set. This dual change of measure can be viewed as 
a dynamic importance sampling algorithm, as introduced formally in Dupuis 
and Wang (2005); i.e., the change of measure depends on the outcome of the 
simulation experiment. 

The motivation for our algorithm is the observation that the SFPE (1.3) 
induces a forward recursive sequence, namely, 

(1.4) Vn = Anmax{Dn,Vn-i)+Bn, n = 1,2, . . . , Vo = v, 

where {{An, Bn, Dn) ■ n G Z+j is an i.i.d. sequence having the same law 
as {A,B,D). Then {Vn} is a Harris recurrent Markov chain, and hence it 
returns with probability one to any set C intersecting the support of its sta- 
tionary measure. Thus we may study {Vn} over the excursions from a set C. 
It is important to observe that in many applications, such as with perpetu- 
ities, the mathematical process under study is usually obtained through the 
backward iterates of the given SFPE (as described Letac (1986) or CoUam- 
ore and Vidyashankar (2011), Section 2.1). In particular, the linear recursion 
/(f) = Av+B induces the backward recursive sequence or perpetuity sequence 

(1.5) z.. = V„+|l + -A_ + ... + -j-^, ,, = 1,2,.^^. 

However, {Zn} is not Markovian, and it is much less natural to simulate the 
sequence {Zn} as compared with the forward sequence {V^}- Thus, a central 
aspect of our approach is the conversion of the given perpetuity sequence, via 
its SFPE, into a forward recursive sequence which we then simulate. We note 
that the ruin problem with investments is similar, where again the underlying 
process is typically described via a backward sequence. 

We remark that another importance sampling algorithm, specifically for 
perpetuities with positive i?i's, was recently proposed in Blanchet et al. 
(2011). Their approach is considerably different from ours, as they analyze 
(1.5) directly, and thus their methods do not generalize to the other processes 
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studied in this paper, such as the ruin problem with investments or further 
processes governed by Letac's Model E. Moreover, in the context of perpetu- 
ity sequences, their methods yield quite different results compared with ours; 
see the discussion in Section 4 below. In another direction, although distinct 
from the main focus of this paper, there is also an interesting connection of 
our work to the problem of simulating the stationary distribution of perpetu- 
ities, which has received some recent attention in the literature (c/. Fill and 
Ruber (2010) or Devroye and Fawzi (2010)). 

We conclude this introduction with a brief discussion of our results and 
some of the main contributions of this article. Section 2 is devoted to a de- 
scription of our main results. Utilizing a dual change of measure over a random 
time interval and the forward recursive sequence {Vn} generated by the given 
SFPE, we obtain an algorithm which, based on our main theorems, is: (i) con- 
sistent; (ii) asymptotically efficient in the sense that it has bounded relative 
error; and (iii) optimal among a wide class of possible importance sampling 
algorithms. 

In Section 3, we provide a proof of consistency and asymptotic efficiency 
for our algorithm. The proof of consistency requires that we relate P {V > u) 
to the number of exceedances of {V^} over a cycle which emanates from a 
given set C C M and terminates upon the return to C. The representation 
formula we obtain involves an embedding of the given Markov chain into a 
Markov additive process and uses aspects of Markov renewal theory. 

In the proof of efficiency, we encounter perpetuity sequences similar to 
(1.5), albeit more general. In the context of (1.5), it is known that while 
the series converges to a finite limit under minimal conditions, the necessary 
and sufficient condition for the L^ convergence of {Zn} is that E [^~ ] < 1; 
cf. Alsmeyer et al. (2009). However, our analysis will involve moments of 
quantities similar to {Zn}, but where E [^"'^l is greater than one, and hence 
our perpetuity sequences will necessarily diverge in L^. To circumvent this 
problem, we will study these perpetuity sequences over randomly stopped 
intervals in the sense described in the previous paragraph, that is, over cycles 
emanating from, and returning to, a given subset C of M. Thus, instead we 
study Zk, where K denotes the first return time to C, and we establish the Lp 
convergence of this quantity. We establish this convergence not only under 
the assumption that E [^4"^^] € (l,oo) but, more surprisingly, also when 
E [^"''^j = oo. This result, obtained in the proof of Theorem 2.3 below, will 
rely on estimates for the return times of {V^} to C. It is worth noting that if K 
were replaced by the more commonly studied regeneration time r of the chain 
{Vn}, then the existing literature on Markov chain theory would not shed 
much light on the tails of r, and hence the convergence of Vr- Thus, the fact 
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that K has sufficient exponential tails for the convergence of Vk is due to the 
recursive structure of the particular class of Markov chains which we consider 
here and seems to be a general property for this class of Markov chains. 
These results concerning the moments of L^-divergent perpetuity sequences 
complement the known literature on perpetuities and appear to be of some 
independent interest. 

In Section 4, which contains examples, numerical results and a discussion 
of the implementation of the algorithm, we describe various numerical strate- 
gies for dealing with models of practical interest, such as ruin models with 
stochastic investments and the ARCH(l) and GARCH(1,1) financial time 
series models. 

Sections 5 and 6 are devoted to sharp asymptotic estimates for the running 
time of the algorithm and to the optimality of the algorithm, respectively. In 
particular, motivated by the Wentzell-Freidlin theory of large deviations, we 
consider other possible measures for generating the process {Ki}, which are 
allowed to depend on the level of the scaled process {logl/^_i/logn} and on 
whether or not {Vn} has exceeded u prior to the present time. We then show 
that our algorithm is, in fact, the only one which attains bounded relative 
error, establishing rather definitively the validity of the algorithm, at least in 
an asymptotic sense as the tail parameter n — )■ oo. In proving these results, 
particularly relating to the running time of the algorithm, we encounter a 
variety of nonstandard issues which we resolve using techniques involving first 
entrance and last exit times of the Markov chain {Vn} generated under various 
measures, and nonlinear renewal theory. Finally, some concluding remarks are 
given in Section 7. 

2. The algorithm and a statement of the main results. 

2.1. Background, hypotheses, and the algorithm. In this section we intro- 
duce our importance sampling algorithm and discuss its theoretical proper- 
ties. We begin by describing the so-called forward and backward iterates of 
an SFPE. In general, an SFPE can be written as a function of the unknown 
random variable V and an environmental random vector Y, specifically, 

(2.1) V = F{V,Y), 

where F : M x M — )■ M is a deterministic measurable function which is 
continuous in its first component. For convenience, we also write Fy{v) = 
F[v,y), where y G M . In particular, Letac's (1986) Model E — which will be 

the main focus of this article — is given by y = Fy{V), where 

FYiv,Y) = Ama.x{D,v) + B and Y = {A,B,D), 
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where typically {A, B, D) are correlated. 

Let V be an element of the range of -F, and let {Yn} be an independent and 

identically distributed (i.i.d.) sequence of random variables such that Yn = Y 
for all n. Then the forward equation generated by the SFPE (2.1) is defined 
by 

(2.2) Vn{v) = FY^oFY^_,o...oFY,{v). n = l,2,..., V^ = v, 
and the backward equation generated by this SFPE is defined by 

(2.3) Zn{v) = FY,oFY^o...oFYSv). n = l,2,..., Zo = v. 

While the forward equation is always Markovian, the backward equation need 
not be Markovian; however, for every v and every n, Vn{v) and Zn{v) are 
identically distributed. This observation is critical since it suggests that — 
regardless of whether the SFPE was originally obtained via forward or back- 
ward iteration — the natural way to analyze the process is through its forward 
iterates. 

Returning to Letac's Model E, the forward equation reduces to 

(2.4) Vn = Anma-K{Dn,Vn-i)+Bn, n = l,2,..., Vo = v, 
which upon iteration yields (with Bq = Vq) that 



(2.5) K = max| ^S, \[ A,,\/ 

i=0 j=i+l k=l 



i=k j=i-\-l j=k 



The regularity properties of {Kt} have been described in CoUamore and 
Vidyashankar (2011). To state the properties relevant for this article, we first 
need to introduce some notations and hypotheses, as follows. Let 

\{a) = E [A'^] and A(q) = log A(q), a G M. 

Let ^ denote the distribution of 1" := {\og A, B,D)^ and let /Xq denote the 
a-shifted distribution with respect to the first variable; that is, 

(2.6) i^a{E) := T7^ / e""d^(x, y, z), Ee B{W^), aeW, 

A(a) Je 

where here and in the following, B{W^) denotes the Borel sets of M^, for any 

d G Z_|_. Let Eq [•] denote expectation with respect to this a-shifted measure. 

For any random variable X, let 2,{X) denote the probability law of X, and 

let supp {X) denote the support of X. Also, write X ~ 2,{X) to denote that 
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X has this probabihty law. Given an i.i.d. sequence {^j}, we wih often write 
X for a "generic" element of this sequence. Finally, for any function /, let 
dom (/) denote the domain of /, and let /' denote its first derivative, f^"^' its 
second derivative, and so on. 

Hypotheses: 

{Hi): The random variable A has an absolutely continuous component 
with respect to Lebesgue measure and satisfies A(^) = for 
some .^ G (0, oo). Moreover, A*-^) is finite on {0,.^}. 

{H2): E[|S|«] < 00 and ^[{A\D\f] < 00. 

{H3): P{A>l,B>0)>OorP{A>l,B>0,D>0)>0. 

The parameter ^ appearing in (Hi) will play an important role in the 
sequel, since we will use the .^-shifted measure defined in (2.6) to develop 
our algorithm. Under the above hypotheses, it is shown in Collamore and 
Vidyashankar (2011), Section 5, that the following path properties hold: 

Lemma 2.1. Assume Letac's Model E, and let {Vn} denote the forward 
recursive sequence corresponding to this SFPE. Assume that {Hi), {H2), and 
{H^) are satisfied. Then: 

(i) {V^} is if -irreducible and geometrically ergodic. Moreover, 

lim P{Vn>u) = P{V>u). 

(ii) Under the measure fi^, 

Vn /^ +00 w.p. 1 as n — 7- 00. 

Thus, the Markov chain {Vn} is transient in the ^^-shifted measure. 

We note that the proof of part (i), providing the uniqueness of the solution 
of the SFPE, is obtained from Letac's principle. The following tail estimate 
is fundamental for the development of the results of this paper. 

Theorem 2.1. Assume Letac's Model E, and suppose that {Hi), {H2), 
and {H3) are satisfied. Then 

(2.7) lim u^P {V >u) = C 

for a finite positive constant C . 

The limit result (2.7) appears in Goldie (1991), while the explicit identifi- 
cation of the constant C is given in Collamore and Vidyashankar (2011). 
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The above theorem shows that the tail of V decays at a polynomial rate 
and hence for large u, P (y > u) constitutes a rare event probability. Thus 
standard Monte Carlo methods break down for large u, in the sense that 
the relative error of the estimator will tend to infinity as the probability in 
question tends to zero (c/. Asmussen and Glynn (2007), Chapter 6). The es- 
sential reason is that the Markov chain will frequently fail to reach the high 
level u. Large deviation theory then suggests that we consider shifted distri- 
butions, and based on known techniques used in the classical ruin problem, 
it is natural to shift using the parameter ^ > satisfying E [^^] = 1. As a 
starting point, we observe that under the measure /i^, the process {Vn} will 
be transient by Lemma 2.1 (ii). Thus under the shifted measure /U^, {Vn} will 
ultimately attain the high level u. 

Hence a reasonable first step is to simulate the process from the distribution 
fi(^. To relate P {V > u) to the paths of {Vn} under /i^-measure, let C := 
[—M, M] for some M > 0, and let vr denote the stationary distribution of 
{Vn}, and define a probability measure 7 on C by setting 

(2.8) 7(^) = ^^y , for aU Borel subsets E of C. 

7r(C) 

Let K := inf {n > 1 : Vn & C} denote the first return time to C. Then we will 
establish below that 

K-l 

(2.9) F{V>u) = E^ [Nu] , where iV„ := ^ l{y„>n}- 

i=0 

This representation formula suggests that we simulate {Vn} over a cycle em- 
anating from the set C and then returning to C. 

For (2.9) to be useful, it is necessary that we simulate the process over a 
cycle that terminates. Since {Vn} is transient in the ^-shifted measure, it is 
natural to simulate in the original measure after the rare event {Vn > u} has 
occurred. Thus, we propose to simulate the process under a dual change of 
measure, namely 

(D) £{logAn,Bn,Dn) = \ ^^ for u = 1 T^, 

^ ^ [ fi tor n > lu, 

where ^^ is defined as in (2.6) and S, is given as in (Hi). Roughly speaking, the 
dual measure shifts the distribution of log An on a path of random duration 
terminating at time Tu, reverting to the original measure thereafter. Let Ej) [•] 
denote expectation with respect to the dual measure (Ti). 

Finally, to relate the simulated quantity in the ^-shifted measure to the 
required probability in the original measure, introduce an appropriate weight- 
ing factor. In the proof of Theorem 2.2 below, we will show that Ej) [£y] = 
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7r(C)E2) [NuC ^'^^«1|x^<^} 1^0 ~ 7]) where 7 is given as in (2.8). Using this 
identity, it is natural to introduce the importance samphng estimator 

(2.10) £u = Nue~^^^^l{T^^K}- 

Then it wiU be shown in the proof of Theorem 2.2 that tt{C)£u is an unbiased 
estimator for P {V > u). However, since in practice the stationary distribu- 
tion TT and hence 7 is seldom known — even if the underlying distribution of 
(A,B,D) is known — we first run multiple realizations of {Vn} according to 
the known fi and thereby estimate 7r(C) and 7. Let 7rfc(C) denote the estimate 
obtained for 7r(C), and let £u,n denote the estimate obtained upon averaging 
the realizations of £u- This yields the estimator Tfk{C)£u,n- The above discus- 
sion can be formalized as an algorithm as follows: 

Algorithm. Step 1. Estimate the stationary measure vr of the process 
{Vn}- Namely, generate i.i.d. realizations {Kt,j : J = 1,...,A;} oi Vn and 
estimate 

1 ^ 

i=i 
Set ^k{E) = -7rfe(£')/7rfc(C), for any Borel set E C C, and set 7a;(C'^) = 0. 

Step 2. Generate a new process having the initial state Vq ~ 7fc- Now, given 
Vn-i, generate Vn using the forward recursion of the SFPE and the shifted 
measure fj,^. li Vn G C, then stop and set f^ = 0. If y„ > u, then go to Step 
3. If not, then repeat Step 2. 

Step 3. Generate the process {Vn} according to its original distribution 
until the first time the process enters the set C. 

Step 4- Compute the importance sampling estimator, namely 

where N^ denotes the number of exceedances of the process {V^} above level 
u prior to time K. We notice here that £u = £u{k)-, since the initial distri- 
bution of the Markov chain depends on k. (We will frequently suppress the 
dependence on k when there is no cause for confusion.) 

Step 5. Repeat Steps 2-4 n times, each time calculating the importance 
sampling estimator Eu, yielding the j'th realization Suj, j = 1, . . . ,n. Then 
compute the average. 



1 " 

n ^-^ 

Finally, estimate V {V > u) by setting this quantity equal to f[i^{C)£u,nik)- 
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It is worth observing that in the special case that D = 1 and B = 0, Letac's 
Model E reduces to multiplicative random walk, and in this case, one can 
always take 7 to be point mass at {1}, at which point the process regenerates. 
In this much-simplified setting, our algorithm reduces to a more standard 
regenerative importance sampling regime, as may be used to evaluate the 
stationary exceedance probabilities in a GI/G/1 queue. We now turn to the 
properties of our algorithm. 

2.2. Consistency and efficiency of the algorithm. Our first result is con- 
cerned with consistency; that is, for any fixed u, the estimator given in (2.10) 
converges to P (y > u) with probability one (w.p. 1) as w — )■ 00. 

Theorem 2.2. Assume that hypotheses (Hi), {H2), and {H3) hold. Then 
for any C intersecting supp(7r) and for any u sufficiently large .such that u ^ C, 
the algorithm is strongly consistent; that is, 

(2.11) lim lim 7rk{C)iun{k) = 'P{V >u) a.s. 

fc^oo n— >-oo ' 

Remark 2.1. If the stationary distribution vr of {Vn} is known on C (as 
would be the case if C = {v} for some point f G M), then it will follow from 
the proof that Tr{C)£u,n is an unbiased estimator for F (V > u). 

We emphasize that the choice of the set C := [— M, M] is not critical for 
obtaining consistency, and C can be taken to be any set having a non-empty 
intersection with the support of the stationary measure tt. However, under 
an appropriate choice of M, we will establish that the estimator is strongly 
efficient in the sense that it has bounded relative error. Before stating this 
next result, we introduce a slight modification of hypothesis {H2). 

Hypothesis: 

[h'^): E[(^~1|5|2)"] < 00 andE[(yl|i:'|2)"] < 00, for some a > C- 

Theorem 2.3. Assume that hypotheses {Hi), {H2), and (H^) hold and 
that X(a) < 00 for some a < —^. Then there exists an M > such that 

(2.12) sup sup u^^Ej, [f2| y^ ^ ^^] < ^_ 

u>Ok£'L+ 

Moreover, even if \{a) = 00 for all a < —^, then (2.12) still holds provided 
that E [{\D\ + A-^\B\)Y] < 00 for all q > 0. 
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Eq. (2.12) implies that our estimator exhibits bounded relative error, that 
is, the ratio Ej)^^ [£^^] / (E^^j [£u\) is uniformly bounded in u and, in partic- 
ular, remains bounded in the limit as w — )■ oo. It is worth noticing that even if 
A(— ^) < oo, the relative error could not be bounded without the presence of 
the term 1 ry^^^^i in the definition of £u- For further discussion of this point 
and its relationship to perpetuity estimates, see the discussion in Section 3 
following the proof of this theorem. 

A good choice of M is critical for the practical usefulness of the algorithm. 
A canonical method for choosing M can be based on the drift condition 
satisfied by {Vn}, generated in its a-shifted measure, with a = and a = —^. 
Indeed, we will use this approach in Lemma 3.3 to provide a formula for 
M. Admittedly, this method may not always yield a good choice for M. 
For instance, a relatively large value for M could lead to an increase in the 
initial rejections, thus leading to an increase in the number of Monte Carlo 
simulations needed for accurate results. Indeed, the choice of M is problem- 
dependent and we outline an alternative numerical method in Section 4 to 
address this issue. 

2.3. Running time of the algorithm. Next we focus on the running time of 
the proposed algorithm, which will clearly depend on the following three quan- 
tities under the dual change of measure: (i) the return time K on {K < oo}; 
(ii) the time r„ needed to reach the level u; and (iii) the time K — Tu to return 
to the set C, where the latter quantity is evaluated on the set {K > Tu}. The 
expected behavior of these three quantities is the subject of our next main 
theorem. First we need to impose a mild regularity condition. 

Hypothesis: 

{H4): F({yi < l|Vo = f) = o{v^'^) as V ^ CO, for some e > 0. 

Theorem 2.4. Assume that hypotheses {Hi)-{H4) hold. Then: 

(i) Under the measure nxi, the return time K has finite expectation on 
{K < oo}; that is, Ej) [KIj^^^qo}] < oo. 

(ii) Conditional on the event {T„ < K}, we have under the measure ^uj) that 

1 



E- 



'D 



logu 



Tn < K 



A'(0 



as ti — 7- oo. 



(iii) Furthermore, conditional on the event {T„ < K}, we also have that 

„ \ K-Tu ^ ^jA ^ 1 ^ 

Ej, — Tu<K ^ TTTTTvu as w ^ oo. 

[ logu J |A'(0)| 



12 J. F. COLLAMORE ET AL. 

Note that the ultimate objective of the algorithm is to minimize cost, i.e., 
the total number of Monte Carlo simulations needed to attain a given accu- 
racy. Thus, it is upon combining Theorem 2.3 with Theorem 2.4 that we see 
that our algorithm is actually efficient. For further details, see Remark 2.2 
below. 

Several sufficient conditions for (-^3) can be provided. For the case B > 
and D = (which corresponds to the special case of a perpetuity sequence), 
it is easy to see that (-^3) always holds, since by Chebyshev's inequality, 

(2.13) P5(Vi < 11^0 = v)<P{A< v-^) < Eg [A-'] v'', 

and for e > sufficiently small. Eg [A'"] = A(^ - e) G (0, 1). When B and 
D assume values in R, a slight modification of (2.13) yields that (-^4) holds 
under appropriate moment conditions on B/A and D. In general, (-^4) is a 
very weak condition which is easily verified in all of the examples considered 
in this article. 

2.4. Optimality of the algorithm. We conclude with a comparison of our 
algorithm to other algorithms obtained using alternative measure transforma- 
tions. A natural alternative would be to simulate with a measure ^a until the 
time Tu = inf{n : Vn > u} and revert to some other measure /i^ thereafter. 
More generally, we may consider a general class of distributions with some 
form of state dependence, as we now describe. 

Let v{-;w,q) denote a probability measure on ;B(M^) indexed by two pa- 
rameters, w G [0, 1] and q € {0, 1}. Intuitively, {w,q) denotes a realization of 
iW^,Qn), where 

, logK-i 
^n-= ^pg^ and Qn:=l{T^<n}- 

Now suppose that we shift the driving sequence {Yn} of the process {Vn}, i.e., 
we change the distribution of Yn := {An, Bn, Dn) ~ /i. In a general sense, we 
would like to allow the simulation distribution of Yn to depend on (VF„, Qn), 
i.e. on the level of the scaled process {log V„/logu} at the prior time n — 1, 
and on whether or not the process {Vn} has exceeded the level u by that time. 
To be more precise, let W^ be defined as above and set 

W =1 ^n, if<G[0,l], 

'^ \ (W^Al)VO otherwise. 
Next, suppose that {I'n} is a family of random measures satisfying 

(2.14) i^niE) = u{E; Wn, Qn), for aU E G B{M.^), 
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for some given measure u. 

Condition (Co): For any measure z/, where fi <^ v, 



Es 



log(^iYn;Wn,Qr. 



Wn = W 



is piecewise continuous as a function of ttJ, for any fixed Qn G {0, 1}. 
Class mi: Define 

(2.15) DJt = < {fn} '■ ^n{-) = ^(s Wn, Qn) and u satisfies (Cq), for all n >. 

Thus, 931 denotes the class of all sequences of measures {vn} obtained from 
the construction (2.14), where (Co) is satisfied and we also assume that ^ 
is absolutely continuous with respect to v[-;w,q) for any w G [0, 1] and q G 
{0, 1}. Set V = {z^i, 1^2, • • •}, and let Ej^ [•] denote expectation with respect to 
this state-dependent collection of simulation distributions. 

From a practical perspective, (Co) is no real restriction on the class of 
possible simulation distributions. In practice, it is natural to simulate using a 
class of exponential transforms determined by a multidimensional parameter 
j3 = f3{w,q), where for some normalizing factor c, 

cJe 

and -C(log An, Bn, Dn) ^ f^/3- In other words, we shift all three members of the 
driving sequence Yn = (log A^, -Bn, Dn) in some way, allowing dependence on 
the history of the process through the parameters {w,q), where w corresponds 
to the realization of log V^_i/logti while q corresponds to the realization of 
'^{Tu<n}- It is reasonable to expect that the optimal P{w, q) depends on w, or 
that it changes in the event that a threshold is achieved, e.g., when the process 
Vn first exceeds u. Such scenarios will always lead to members v = {vi, 1/2, . . ■) 
belonging to the class Tl and, in particular, satisfying (Co). 

Now suppose that simulation is performed using a modification of our pre- 
vious algorithm, where Y^ ~ f„ for some v = {i^i, zv2, • • •} G 93T. For the 
importance sampling estimator, set 

(2.16) 4-) = ivJ n ^(>1; W„ qA 1{t.<k}. 

Let TTfe denote an empirical estimate for tt, as obtained in Step 1 of our main 
algorithm, and let S^^l , . . . , Su'^n denote simulated estimates for £u obtained 
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by repeating this algorithm but with {^'n} in place of the dual measure (2D). 
Then it is easily shown by a slight modification of Theorem 2.2 that 

(2.17) hm ^fc(„)^a = P (V^ > ^) , 

where S^^n denotes the average of n simulated samples of £u ■ It remains to 
compare the variance of these estimators, which is the subject of the next 
theorem. 

Theorem 2.5. Assume that the conditions of Theorems 2.3 and 2.4 hold. 
Let u be a probability measure on B(M^) indexed by parameters w £ [0, 1] and 
q £ {0, 1}, and assume that u G 9Jt. Then for any for any initial state v £ C, 



(2.18) liminf -^ log (u^^E^ \ [Sl^^] 

u^oo log U \ L ^ ' 



Vb 



> 0. 



Moreover, equality holds in this inequality if and only if iy(-;w,0) = /i^ and 
v{-] w, 1) = n, for all w G [0, 1]. Thus, the dual measure in (D) is the unique 
optimal simulation strategy in the class Tl. 

Remark 2.2. In practice, the objective should be to minimize the total 
number of random variables generated in order to obtain a given accuracy. 
This grows according to 



(2.19) Var(4'^)) {ciE, [K |T„ <K] + C2E, [^l{T„>i^}] } as 



00 



for appropriate constants ci and C2; cf Hammersley and Handscomb (1964), 
Siegmund (1976). However, as a consequence of Theorem 2.4 (ii) and (iii), we 
have under the dual measure (2)) that 

E-x) [K \Tu < K] ^ Qlogu as li — ^ 00, 

for some positive constant Q, while the second term in (2.19) converges to 
a finite constant, by Theorem 2.4 (i). Thus, under the dual measure, it is 
sufficient to minimize the second moment in the asymptotic limit as u — )■ 00. 
Consequently we conclude that simulation under the dual measure (S) is 
indeed asymptotically efficient and optimal. 

3. Proofs of consistency and efficiency. Set Aq = 1, and define 

n 

Sn = Y,logAi, i = l,2,...; So = 0; 



1=1 

V - °° B 

Zn = —. . , n = 0,l,...; and Z^^' =y — . l/R-^ni; 
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where 

Bo = \Vo\ and Bn = An\Dn\ + \Bn\- 

In the following lemma, we summarize some regularity properties satisfied 
by {Zn} and the corresponding perpetuity sequence Z^P'. 

Lemma 3.1. Assume Letac's Model E, and suppose that {Hi), {H2), and 
(H^) are satisfied. Then: 

(i) Zn — > Z with respect to the measure fi^, where Z is a proper random 
variable supported on (0,cxd). 

(ii) E^[(Z(P))^] < 00, and for all n and all u, 

(3.1) \Znl{n<K}\<Z^''^ and \Zt^1{t^<k}\ < Z^'^ ■ 

Moreover, 

00 ^ 

(3-2) ZT„l{T„<i^} < ^ -. ^^-l{n<T„<ii-}- 

^0 ■ ■ ■ ^n 
n=0 

Proof. These properties through (3.1) are established in CoUamore and 
Vidyashankar (2011), Lemma 5.5. (Although we have replaced the stopping 
time r of that article with K in (3.1), the proofs remain unchanged.) 

For (3.2), apply Eq. (5.32) in Lemma 5.5 of Collamore and Vidyashankar 
(2011) to obtain that 

n ~ 

(3.3) ^n < E 4 \ ' n = 0,l,.... 

Then (3.2) follows from (3.3). D 

A second basic result which will be used throughout this article is a change 
of measure formula, which allows one to compare an expectation under a par- 
ticular measure to an expectation under another measure. Frequently we will 
take this other measure to be the dual measure. In the following lemma as- 
sume, quite generally, that Un denotes the distribution of {(y4„, Bn, -D„)}, pos- 
sibly depending on the past history of {Fn}, that is on the vector (Vq, • • • , V^n-i)- 
To be more precise, let v„ = (uq, . . . , Vn-i) denote a vector in M", and intro- 
duce the family of measures 

{vn{E■^rn) : E G S(M^), v„ G M",n G Z+}. 

Then our objective is to compare expectations under the original measure to 
expectations under the driving sequence 

y„ := {An, Bn, Dn) ~ f„(-; Vq, . . . , K-l). 
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Let E^(„) [•] denote expectation with respect to the family {ui, . . . , Vn} just 
described. With a sUght abuse of notation, we will surpress the dependence 
on n in this last expectation. Then we have the following: 

Lemma 3.2. Assume Letac's Model E, and suppose that ^ <^ Vn a.s. 
for all n. Let g : M°° — )■ [0, oo] he a deterministic function, and let gn de- 
note its projection onto the first n + 1 coordinates; that is, gn{xo, . . . , x„) = 
g{xo, • . . , Xn, 0,0, . . .). Then for any n £ N, 



(3.4) B[gniVo,...,Vn)] = B^ 



\i=l 



YlT^(^^n9n{Vo,...,Vn 



Proof. This result follows by induction. Assume that the result holds for 
some positive integer n, for any given function g, and consider 



(3.5) 



E. 






E,, 



fl^{Y,)]hn{Vo,...,Vn 



\i=l 



dvi 



where 

/i„(yo,...,K):=E, 



dfi 

dVn+ 



■(^n+l)5'n+l(Vb, • • • , Vn+l) 



1^0, • . . , K 



By direct calculation, 

/i„(Vb, . . . , K) = E [gn+i{Vo, . . . , K+i)| Vo, • • • , K 
Now apply the inductive hypothesis to obtain 



E. 



d/U 



\{-r{Y^)]hn{Vo,■■■,Vn) 



\i=l 



E 



/i„(Vb,...,K) 



Substituting these last two equations into the right-hand side of (3.5) yields 
(3.4), as required. □ 

We will be particularly interested in applying the previous lemma to dual 
measures of the form 



i^c 



^{\og An, Bn,Dr, 



fia for n= 1,. ..,Tu, 
fi for n > Tu, 



where a G dom (A). By conditioning on {T^ = m, K = n} and summing over 
all possible values of m and n, we obtain the following: 
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Corollary 3.1. Assume the conditions of the previous lemma. Then for 
any a G dom (A), 

(3.6) E [gKiVo, ...,Vk)]= Ej,„ [<7x(Fo, • • • , VK)e-°^^- (A(a))^" 1{T„<X}" 

+ Ej,^ [gKiVo, ..., yx)e-"^^ (A(a))^ 1{T„>X}" , 

where E'X)^[-] denotes expectation with respect to the dual measure in (Sq-) 
above. 

For notational convenience, here and in the fonowing, we write Ej) [•] as 
Ej) [•], i.e., we suppress ^ in the special case that a = $,. Moreover, we use 
the notation Pq, (•) and E^ [•] when the random variables under consideration 
have the distribution Ha for all n. 

We now turn to the proof of Theorem 2.2. The main complications in the 
proof are that: (i) we work with the return times rather than the regeneration 
times of the Markov chain; and (ii) we utilize an empirical approximation to 
the stationary distribution. The consequence of (i) is that the cycles induced 
by the return times are not i.i.d. We begin with a summary of the main steps 
of the proof. 

Step 1. Letting Ki, K2, . . . denote the return times of the process {F„} to 
C, we show that the stopped process Xn = Vk^, n = 1, 2, . . . , is a stationary 
Markov chain whose stationary distribution is given by ^{E) = it{E)/it{C), 
for any Borel set E C C, where vr is the stationary distribution of {y„}. 

Step 2. We derive an expression for the total number of exceedances N^ 
above the level u which occur over a cycle emanating from C and terminating 
upon its return to C, assuming that the initial distribution satisfies Vq ~ 7. 

Step 3. We relate the expression obtained in Step 2 to the object of interest, 
namely P {V > u). 

Step 4- We verify that, if the true stationary distribution 7 is replaced with 
the empirical distribution ^1., then the resulting estimator is consistent in the 
asymptotic limit where 7a; =^ 7. 

Proof of Theorem 2.2. Step 1. Let Ki,K2,... denote the successive 
return times of {Vn} to C; that is, Kq = and Kn = inf {i > Kn-i : Vi £ C} 
for n = 1, 2, . . . . Set 

Xn = VK„, n = 0, 1,.... 

Note that {^n} is a Markov chain with state space C. We will now show 
that {X„} has a stationary distribution given by jiE) = ■k{E)/ti{C), for all 
Borel sets E C C. First observe that {V^} is 99- irreducible and geometrically 
ergodic, by Lemma 2.1 (i). It follows by the definitions of irreducibility and 
ergodicity that {X„} is also 99-irreducible and geometrically ergodic. Thus 
{Vn} and {Xn} are both positive chains. 
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Hence applying the law of large numbers (Meyn and Tweedie (1993), The- 
orem 17.1.7 with f{x) = 1e{x), where 1e is the indicator function on the set 
E), we obtain that 

1 " 

(3.7) lim - V IrygEi = 7r(E) a.s., for ah EcC. 

n— >-oo n ^ — ' * 

Next, let 0T„ := Y17=i '^{Xt^c} denote the number of visits of {X„} to C occur- 
ring by time n. Since {Vn} is a recurrent Markov chain, it follows that 9T„ f oo 
w.p. 1. Applying (3.7) with E = C, we obtain that tt(C) = lim„^oo ^n/n a-s. 
Similarly, applying (3.7) with E C C and utilizing the definition of {Xn}, we 
also obtain that 

Tr{E) = lim -^ I -^ > ^ l|x.6£| 1 a.s. 



Hence 



n— >oo fi 



m„ 



[^^EMx.mj 



i=l 

Finally, by yet another application of the law of large numbers for Markov 
chains, now applied to {Xn}, we obtain that the left-hand side of (3.8) con- 
verges to the stationary distribution of {Xn}- Thus, it follows that jiE) := 
tt{E)/tt{C) is the stationary distribution of {Xn}- 

Step 2. Let Nu^i, Nu^2, ■ ■ - denote the number of exceedances above level 
u which occur over the successive cycles starting from C; that is, N^^n '■= 
E^xii Mv.>u}- Set S^ = and 5^ = 7V„,i + • • • + iV„,„, n = 1, 2, . . . . 
Then {{Xn, Sn) : n = 0, 1, . . .} is a Markov random walk. To view {Sn} as a 
sum of the functionals of a Markov chain, first introduce the adjoined Markov 
chain {{Xn, Nu,n)}- Since the distribution of Nu^n is determined entirely by 
Xn, where {Xn} has stationary distribution 7, it follows that the stationary 
distribution of this adjoined chain is given by 



(3.9) 7(^ X F) 



[ P (iV„,i G F|Xi = x) d-f{x) 
Je 



Moreover, as {-^n} is a positive Harris chain, it follows that {[Xn, Nu,n)} 
is also a positive Harris chain. Hence the law of large numbers for Markov 
chains applies, yielding (after observing E^ [Xu] = E^ [Xu]) that 

gN -^ K„-l 

(3.10) E^[iV„] = lim -^ ■-= lim - V l|y^>„} a.s., 

n=0 
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provided that the expectation on the left-hand side is finite. 

To see that this quantity is finite, recaU that {Vn} is geometrically recur- 
rent. Then since {Kn} denotes the return times of this process to its C-set, 
we have that sup^g^ E^, [Kn — Kn-i] < oo. Since Nu,n denotes the number 
of exceedances above level u occurring over the nth such cycle, we obviously 
have Nu^n < (Kn — Kn~i), and thus E^[A'"ul < oo. 

Step 3. Next, we identify the limit on the right-hand side of (3.10). First 
recall that 9T„ is the number of returns to C occurring by time n; thus K'yx^ is 
the last return time occurring in [0,n]. Now apply the law of large numbers 
for Markov chains to obtain that 

(3.11) P(y>^):=7r(^,oo) = Jim -^ J] l^W + E V.>«} • 

We claim that the last term on the right-hand side is asymptotically neg- 
ligible. To justify this fact, we apply a Markovian renewal theorem given in 
Iscoe et al. (1985). To this end, first introduce the augmented chain {iVn-, rjn)}, 
where {r/„} is an i.i.d. sequence of Bernoulli random variables, independent of 
{!/„}, with P (r/„ = 1) = (5. Here 5 is the constant appearing in the minoriza- 
tion condition of the Markov chain {T^}. Then the event {Vn G C,?7„ = 1} is 
a regeneration time (Nummelin (1984), Chapter 4). Let I{n) denote the last 
regeneration time occurring in the interval [0, n], and let J{n) denote the first 
regeneration time occurring after this time. Since regeneration only occurs at 
the return times to C for which rji = 1 (so regeneration is "less frequent" than 
returns to C), we have that /(n) < Ktji^ < n < J{n) — 1. Hence 

n J(n)-1 

i=K<yin i=I(n) 

Now apply Lemma 6.2 of Iscoe et al. (1985) to the right-hand side of this 
equation. Letting r denote a typical interregeneration time and letting zv 
denote the measure appearing in the minorization for {V^}, then Lemma 6.2 
of Iscoe et al. (1985) states that 

J(n)-1 ^ 

(3.13) lim Y lW>«} = iMl^ 



n- 

i=I{n) 



i=0 



Vn ^ f 



Since {Vn} is geometrically ergodic, E [r] < oo. Substituting (3.12) and (3.13) 
into (3.11), we conclude 



(3.14) 
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Since 9T„/n — )■ it{C) as n — t- oo, it follows from (3.10) and (3.14) that 

(3.15) P{V>u) = 7:{C)E^[Nu]. 

Next apply Corollary 3.1 to the second term on the right-hand side of 

(3.15) to obtain that 

(3.16) F{V>u) = 7r(C)Ej,,^ [A^„e-«^^" ^{t^<k}\ ■= vr(C)Es,,^ [£u] ■ 

This shows that if we could choose Vq ~ 7, that is, if we knew the stationary 
distribution of {Vn} exactly, then -k{C)£u would be an unbiased estimator for 
P{V >u). 

Step 4- Now suppose that we approximate 7 with the empirical measure 
7fc described in Step 1 of the algorithm. Then it remains to show that 



(3.17) 


lim E2) 

A; -^00 


N 


ue-«^^"l 


[Tu<K} 


Vo ~ 7fc 








= Ex, 


Nue-^^^-l{T^^K} 


Vo ^ J 


Set 




(3.18) 


H{v) = E 


'D 


Es, [Nu 1 


drje- 


^^-1|T„ 


<K} 


Vo = v 



We now claim that H(v) is uniformly bounded in f G C. To establish this 
claim, first apply Theorem 4.2 of Collamore and Vidyashankar (2011) to ob- 
tain that 

(3.19) E^[Nu\dTjl{T^<K} < (Cl(n)log(^) +C2(U)) l{T^^r}, 



where t > K is the first regeneration time, and Ci{u) 
Moreover, by the definition of {Zn}, 



Cj < 00 for i = 1, 2. 



(3.20) 



-i^Tu 



Vt 



< 



^T„- 



Substituting the last two equations into (3.18) yields 



(3.21) 



\H{v)\ < BEj, 



\yi -\ I 



Vn = V 



for some finite constant Q, independent of u. Next apply Lemma 5.5 (ii) of 
Collamore and Vidyashankar (2011) and the definition of Z^P' to obtain that 



(3.22) 



E 



D 



l7« 1 I 



Vo = v 



<Es, 



^z(p)f 



Vn = v 
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By Lemma 3.1 (ii), the right-hand side of (3.22) is bounded, uniformly in 
w G C. Consequently H{v) is bounded in v £ C. Moreover, the function 
H is also continuous, as can be seen by observing that for any fixed /, 
{Vi{'w), . . . , Vi{w)} converges pointwise to {Vi{v), . . . , Vi{v)} as the initial 
state tf — 7- f . Summing over all realizations {K = 1} and using the definition 
H{v) = Ej) [iVue^^'^^"l|y^</^||Vo = v] obtains continuity. Finally, since 7,^ 
and 7 are both supported on C, it then follows from the weak convergence 
Ik ^ 1 that 

lim [ H{v)d-fk{v) = [ H{v)d-f{v), 
k-^^ Jc Jc 

which is (3.17). D 

Next we turn to the proof of efficiency. Suppose now that Yn ^ ^p for all n, 
where Y = {An, Bn, Dn) and j3 G dom (A), and where ^p is defined according 
to (2.6). Let 



A^(a) = / e°'^d^p{x,y,z), A^(q;) = log A/3(a), for all a G M, 

and note by the definition of ^p that 

(3.23) A^(q) =A(a + /?)-A(/3). 

In the following lemma, we summarize some standard results concerning 
the return times to the set C under this /3-shifted measure. Letting P denote 
the transition kernel of {Kt}) then we say that {Vn} satisfies a drift condition 
if there exists a function /i : M — )■ [0, 00) such that 

{V) h{y)P{x,dy)<ph{x), for all x ^ C, 

where p £ (0, 1) and C is some Borel subset of M. We note that this definition 
differs slightly from the more standard definition given in Meyn and Tweedie 
(1993), but will be more convenient for our purposes here. 

Lemma 3.3. Assume Letac's Model E, and suppose that {Hi), {H2), and 
{H^) are satisfied. Let {Vn} denote the forward recursive sequence generated 
by this SFPE under the measure pp, chosen such that infQ,>o Xi3{a) < 1. Then 
the drift condition (D) holds with h{x) = |x|", where a > Q is any constant 
satisfying the equation A^(a) < 0. Moreover, in this case, we may take p = pp 
and C = [—Mp,Mp\, where 



(3.24) pp := tXp{a), for some t G f 1 



Xp{a) 
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and 



(3.25) Mp := < 
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B^j[B"]j (^A^(a)(t-l)j , if a G (0,1), 

^(e^S-])'^" ((A^(a))i/° (tV" -i)y\ if „ > 1. 
Furthermore, for any {pp^Mp) satisfying this pair of equations, we have 



(3.26) 



supP/3 {K > n\VQ = v) < p'p, for all n £ Z+. 



Proof. To verify (T>) and identify p^ and Mjj explicitly, first observe that 

(3.27) \Vn\ < An\Vn-i\ + B^, for all n, 

where Bn '■= An\Dn\ + \Bn\- If a > 1, then it follows by Minkowskii's inequal- 
ity that 

(3.28) Ep[\Vir\Vo = v] < ((E4^°])^/%+(E45"])'^"y 

'Mb 



Pp^ 



tV" 



+ 



l/a^ 



1/a 
Pp ^ 



where 



tXfiia). 



Thus (V) holds with h{v) = |t;|" and p/B = t\p{a) < 1 for i e (l, (A/3(a))"^), 
provided that the quantity inside the parentheses on the right of (3.28) is less 
than one for v ^C. Setting t~^/" + (E^fS"]") /{pgv) = 1 and solving for 
V, we obtain the expression given on the right-hand side of (3.25) for the case 
a> 1. 

If a < 1, then an analogous condition is obtained by employing the deter- 
ministic inequality |x + y|" < |x|" + |y|° in place of Holder's inequality. 

Once {T>) is established, the proof of (3.26) is obtained as in Nummelin 
and Tuominen (1982) or Nummelin (1984). Namely, observe by (P) and an 
inductive argument that 



(3.29) 



E 



^(K) n '^{v^ic} 



Vn 



1=1 



<p^/i(yo), n = l,2, 



Since h{Vn) > h{Vo) on {Vn (^ C, Vq £ C}, it follows that P {K > n) < pj^, 
uniformly in Vq G C. □ 

Proof of Theorem 2.3. In the proof, set 



Pn = An 



■A 



n = 0, 1,..., 
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and assume Vq = v £ C. We will show that the result holds uniformly in 
veC. 

Case 1: A(a) < oo, for some a < —^. To evaluate 

Ex, [£l] :=Ex, [iV2e-2«^-"l|r„<;,| 

note by definition that Vne~^" := Vn/Pn '■= -^n- Since Vr^ > u, it follows that 
< ue~^'^^ < Zt^. Consequently by Lemma 3.1 (ii), 



(3.30) 



?«Es, [£l] < Ex, 



^'u E 



oo ^ \ 25 
On \ 



\n=0 



J~^n, 



^{n<Tu<K} 



where {Bn} is defined as in the discussion prior to Lemma 3.1. If 2^ > 1, then 
using Minkowskii's inequality on the right-hand side, we obtain 



(3.31) 



n=0 V 



iV,^(I^V'l 



p / -{n<T^<K} 



1/2? 



E(e 

n=0 



^IPfi ^Bl^'^{n<Tu<K} 



1/2? 



where in the last step we have used the change-of-measure formula in Corol- 
lary 3.1. Now apply Holder's inequality to the right-hand side. Using the 
independence of (yl„,i?„) from l{n_i<T„Ai^}) ^^ obtain that the left-hand 
side of (3.31) is bounded above by 



Y.i^[m)' '[^ 



n=0 



n n 



s( 



l/2s5 



E 



Pn-l'^{n-l<TuAK} 



l/2si 



where r~^ + s^^ = 1. Now set C = s£, for the remainder of the proof. Next, 
by Lemma 3.2 (c/. Corollary 3.1), the last term on the right-hand side of the 
previous equation may be expressed as 



E 



-C 



-l^{n-l<TuAK} 



,n-l 



(A(-C))" 'P-cin-KT^AK). 



On the right-hand side, we have used the fact that for all n < T^, the dual 
measure (Xi_(^) agrees with the (— C)-shifted measure. Substituting the last 
equation into the upper bound for (3.31), we conclude that 



1/2? 



(3.32) {u^^E^[£l]) <Y,Jn{{X{-OyP-c{n-l<TuAK)} 



1/2C 



n=0 
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where 



2r-\\l/'2r( /„[/._! ^ox^lX^/^C 



Jn:={B[N'/]y'^'^(B[{A;;'Blf\) \ n = 0,l,.... 

Since N^ < K, Lemma 3.3 (with (3 = 0) yields 

(3.33) supE [NJ^' \Vo = v~\ < oo, for any finite constant r. 

Moreover, for sufficiently small s > 1 and C = sS,, it follows by {H2) that 



E 



-1r2\C 



{A-^B^) 



< 00. Thus, to show that the quantity on the left-hand side 
of (3.32) is finite, it is sufficient to show that for some C > '^ ^^^^ some t > 1, 

(3.34) P_(^ (n - 1< T„ A K) < (tA(-C))""^\ for all n > some Nq, 

and that this last equation holds uniformly in u and uniformly in v £ C. Note 
that {Tu AK > n- 1} C {K > n - 1}, and by Lemma 3.3, 

(3.35) supP_^ {K>n-l\Vo = v)< (a(-C))""^\ 
vec 

where C := [—M,M] and M > M_^. (Since C > ■? was arbitrary, we have 
replaced M_^ with M_^ in this last expression. We note that we also require 
M > Mq, so that (3.33) will hold.) This proves (3.34), thus establishing the 
theorem for the case 2,^ > 1. 

If 2^ < 1 then, in place of Minkowskii's inequality, we use the determin- 
istic inequality |x -|- t/l" < Ixj" -|- |y|", for a G (0,1]. This yields (3.31), but 
without the powers "1/2^" on the left- and right-hand sides. Once this mod- 
ification has been made, the previous argument can be repeated to obtain 
that u^^Ej) [iS^] is bounded uniformly in u and v £ C, completing the proof. 

Case 2: A(-C) = 00 for C > ^, but E[(^"ii^)"] < 00 for aU a > 0. 

We modify the previous argument, now employing a truncation argument 
as follows. First assume 2^ > 1. Then as before, we obtain that (u^^Ej) [iS^] ) 
is bounded above by the right-hand side of (3.31), and it is sufficient to show 



thatE 



Pn-l'^{n-l<TuAK} 



< 00 for some C > 6 Set Wn = Pn-A 



l^{n-l<TuAK}i 



and first observe that E [Wn] < 00. Indeed, using the inequality (3.27), namely 
\Vn\ < An\Vn-i\ I 1+ , ^; I |, n = l,2,..., 

\ An\Vn-l\ J 

and observing that n — 1 < T„ A i^ =^ \Vi\ G (M, u) for i = 1, . . . , n — 1, we 
obtain as a rough upper bound that 

(3.36) A~^ kI^X [l + ^^\ , i = l,...,n-l on{n-l<r„ Ai^}. 




IMPORTANCE SAMPLING FOR SFPE 



25 



Then this equation yields an upper bound for Pn~i- Using the assumption 
that E {A-^BY < oo for all a > 0, we conclude by (3.36) that E [W„] < oo. 
Next let {Lk} be a sequence of positive real numbers such that L^ | as 
fc — >■ oo, and set Ff. = 0^=1 {^i — ^k}- Assume that L^ has been chosen 
sufficiently small such that 



(3.37) E [WklF^] < 

Then it remains to show that 



1 



A; = 1,2, 



(3.35 



5^E[WfclFj<00. 



fc=0 



To verify that (3.38) holds, set Aq^/. = 1 and 

^n,k = -4„l{A„>Lfc} + ^fcl{A„<Lfe}) 



n 



1,2, 



Let Afe(a) = E[ylf;,], and define Wk = {Ao---Ak-i) ^^{k-KTuAK}- Since 
{k — 1 < Tu A K} C {A; — 1 < K}, it follows by Lemma 3.2 (similarly to 
Corollary 3.1) that 



(3.39) 



E [Wk] < (Afc(-C))'=-^E_c [l{K>k-i}lF,] . 



To evaluate the expectation on the right-hand side, we apply a direct argu- 
ment to obtain a drift condition (rather than Lemma 3.3, as we used before). 
In particular, start with the forward recursion 



(3.40) 



''rij/cl _ -^'n.-klvr 



n,k\ *^n—l,k\ 



1 + 



-Bn 



^n,k\ •'n— l,fc 



n 



1,2,.... 



Write E_^^^ [•] = E_^ [-l^o./fc = w]. Now let /3 > and take the /3th moment 
in (3.40) for the one-step transition which starts at Vq^a; = w. To simplify the 
expression on the right-hand side, use the definition (2.6) of //_^ (but applied 
to the truncated random variable Ai^/.) to express this expectation in terms 
of the original measure. This yields 



(3.41) 



E_ 



■(,w 



\Vi,k\ 



< 



w 



Afc(-C) 



E 



{Ai,k 



< 



1 + 



B, 



Hence by Holder's inequality. 



(3.42) 



E_ 



■(^,UI 



w, 



l,fcl 



< pkw^ t-«E 



1 + 



Bi 

wA 



wAi^k 



9/31 



l/<? 



l,fc 
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where pk '■= [E 



{A 



i,fc) 



MI3-0 



i/p 



t/Xk{-0) and p-i + g~i = 1. 

Set [3 = argmin^ A(q) and choose /3 such that p{l3 — Q) = (3, and assume 
that p has been chosen sufficiently close to one so that the expectation in the 
definition of pk is finite. Since L^ | as A; — )■ oo, 



(3.43) 



lim Afe(a) = A(a), a > 0, 



by a monotone convergence argument. Moreover, {Hi) implies that A achieves 
its minimum on the positive axis, at which point A(/S) < 1. Thus we conclude 
that for t G ( 1, (A(/3)) j and for some constant p G (0, 1), 



(3.44) 



lim Afe(-C)/>fc = t A(/3) 



fc— >oo 



1/p 



<p. 



Then \k{—C,)pk < P for all k > ko, and with this value of p, we obtain by 
(3.42) that 



(3.45) E_<^,^ [\Vi,k\ 
provided that 

(3.46) r'?E 



< 



pw'~ 



Afc(-C)' 



for all k > ko, 



1 + 



B, 



wA 



l,k 



qP 



< 1. 



Our next objective is to find a set C = [—M,M] such that for all w ^ C, 
(3.46) holds. First assume (7/3 > 1 and apply Minkowskii's inequality to the 
left-hand side of (3.46). Then set this quantity equal to one, solve for w, and 
set w = Mfc. After some algebra, this yields 



(3.47) 



Mk 



tV/3 - 1 



E 



5i 



A 



l,k 



qP' 



1/9/9 



-16^9/3 



{A-^B 



clS K 



oo. Using the 



The quantity in parentheses tends to E 

assumption that E[(^~"^i?) ] < co for all a > 0, we conclude that M := 
sup^ Mfc < oo. 

If g/3 < 1, then a similar expression is obtained for M by using the deter- 
ministic inequality |x -|- y\^ < \x\^ + \y\^ in place of Minkowskii's inequality. 

To complete the proof observe that, as in the proof of Lemma 3.3 (c/. 
(3.29)), we obtain upon iterating (3.45) with C = [-M,M] that 



(3.45 



E-c [l{x>fe-i}lFfe] < 



A(-C) 



-k+l 



for all k > ko- 
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Note that on the set F^, the sequences {Vn^k '■ 1 < n < k} and {Vn : 1 < n < 
k} agree, and thus the event {K > k — 1} coincides for these two sequences. 
Finally, substituting (3.48) into (3.39) yields (3.38), as required. If 2,^ < 1, 
then the above argument may be modified in the same manner as for Case 1, 
and the required estimate still holds. □ 

Returning to the remark concerning the assumption that A (a) < oo for 

P„ Bn I = OO (see Alsmeyer et al. 



some a < — ^, we notice that Ej) 
(2009)). Thus, even though the expectation of the perpetuity sequence is 
infinite, it converges when studied on the set {T„ < K}. We emphasize here 
that some of the above calculations could be made much simpler if we were 
to work with assumptions such as i? > and D > 0. However, these are quite 
strong assumptions from an applied perspective. 

4. Examples and simulations. In this section we provide several ex- 
amples which illustrate how to implement the algorithm. As seen in Sections 2 
and 3, while the use of a drift condition provides a formula for M, it may not 
be optimal in a practical sense. This is due to the fact that the estimate for 
V^ typically uses Minkowskii- or Holder-type inequalities, which are usually 
not very sharp. We begin by outlining an alternative method for obtaining 
M and use it to verify that it yields meaningful answers from a practical 
perspective. 

4.1. Numerical procedure for calculating M . The numerical procedure in- 
volves a Monte Carlo method for calculating the conditional expectation ap- 
pearing in the drift condition, that is, for evaluating 



E, 



V(, = v 



E, 



, D \ BY 
Amax ( — ,1 j H j 



when /3 = and /3 = —^. The goal is to find an a such that M := max(Mo, M_^) 
is minimized, where M« satisfies 



E, 



, D \ BY 
Amax ( — ,1 j H j 



< pi3, for all V > M^ and some pjs £ (0, 1). 



In this expression, a is necessarily chosen such that E^ [^"] G (0, 1), and 
hence we expect that p^ G (E^ [^"] ,1). Note that M depends on the choice 
of a; thus, we also minimize over all possible a such that E^ [^°] G (0, 1). 

Let {{Ai, Bi, Di) : 1 < i < N} denote a collection of i.i.d random variables 
having the same distribution as that of {A,B,D). The numerical method for 
finding an optimal choice of M proceeds as follows: 
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Step 1. Using a root finding algorithm, obtain ^. That is, solve for ^ in 
the equation E[^^] = 1. If an analytic expression for E[^^] is not available, 
numerical approximations such as Gauss-Hermite quadrature methods can be 
applied. 

Step 2. For E^ [^"] < 1, use a Monte Carlo procedure to compute E/j [|Vi|"| 
^0 = ^1 and solve for v in the formula 



1 



N 



i=l 



Ai max ( — , 1 H 

V V 



Pp, 



where this quantity is computed in the /3-shifted measure for /? = and 
(3 = —^ and where /j^ < 1. Then select a so that it provides the smallest 
possible value of v. Set Mjs > v fov [3 = and f3 = —^. 
Step 3. Set M = max(Mo,M_^). 

We now turn to some specific examples. 

4.2. The ruin problem with stochastic investments. Assume that the fluc- 
tuations in the insurance business are governed by the classical Cramer- 
Lundberg model, 

Nt 
(4.1) Xt = U + ct-Y,Cn 



n=l 



where u denotes the company's initial capital, c its premiums income rate, 
{Cn} its claims losses, and Nt the number of Poisson claims occurring in the 
interval [0,t]. It is assumed that {Cn} is i.i.d. and independent of {Nt}. We 
now depart from this classical model by assuming that at the discrete times 
n = 1,2,..., the surplus capital is invested, earning stochastic returns {Rn} 
which are assumed to be i.i.d. Let -L„ := — {Xn — Xn-i) denote the losses 
incurred by the insurance business during the nth discrete time interval. Then 
the total capital of the insurance company at time n is described by the 
recursive sequence of equations 

(4.2) Yn = RnYn-i - K, n = 1,2, . . . , Yo = u. 

It is typically assumed that E [log R] > and E [L] < 0. 

Our objective is to determine the probability of ruin using importance 
sampling, namely to estimate 

(4.3) ^(u) := P (y„ < 0, for some n G Z+ \Yo = u) . 
Iterating (4.2) yields y„ = {R1R2 ■ ■ ■ Rn) {Yo — ^n) , where 

(4.4) C^:=^ + ---+ ^" 



Ri RiR2---R, 



■n 
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and hence ^(u) = P (£„ > u, some n). Setting C = (sup^g^ £„) V 0, then 
from (4.4) and an elementary argument we obtain that C satisfies the SFPE 

(4.5) C = {AC + B)+, where A := — and B = — . 

Ri R\ 

This can be viewed as a special case of Letac's Model E with D := —B/A. 

4.2.1. Implementing the algorithm,. To implement the algorithm in this 
example, we generated investment returns according to the Black-Scholes 
model. Specifically, we chose 

(4.6) An = exp \~ \l^~ ^] ~ ^^n \ , for all n, 

where {Zn} is an i.i.d. sequence of standard Gaussian random variables. Then 

(4.7) Aia) = -a(iJ.-^) + 



a^\ a^a^ 



Thus C = 2///0-2 - 1 and ^^ ~ Normal(// - a^/2, a^). 

We set ^ = 0.2 and a^ = 0.25. Regarding the insurance model, we set the 
premiums rate c = 1, {Cn} to be exponential with parameter 1, and {Nt} to be 
a Poisson process with parameter 1/2. Applying the procedure described at 
the beginning of this section to this model, we obtained that Mq = = M_^. 
Thus M = 0. 

In this example, we can actually deduce that M = by a more elementary 
argument. Arguing as in the proof of Lemma 3.3, we obtain that M^ = 

minj=i^2 -^«* ) where 

(4.8) M« = inf , "51'^'".,, , mP = inf "^f'''^;" , 

P "G(o,i)n$ (l-||^i||^_^)V" P ae[i,oo)n$ l-||yli||/3,a 

and $ = {a G M : E^ [A"] < 1}. Here, || • ||/3^q, denotes the La norm under 
the /^-shifted measure up, where we again consider the two cases /3 = and 
/3 = — ^. For each of these cases, this infimum may be computed numerically, 
yielding Mq = = M_^ and thus M = 0, just as before. 

In choosing M = 0, the algorithm in Section 2 simplifies considerably, 
since in this case we can take the measure 7 to be a point mass at the origin. 
Moreover, it can be observed that {0} is an atom of the the Markov chain 
{Vn}- Consequently, it follows that a cycle originating at {0} and then stopped 
upon its return to {0} forms a regeneration cycle of the Markov chain (c/. 
Nummelin (1984), Chapter 4). 
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We implemented the algorithm in Section 2 to estimate the probabilities 
of ruin for u = 10, 20, 100, 500, 10^, 5 x 10^, 10^, 5 x 10"^, and 10^ In 
all our simulations, the distribution in Step 1 was based on k = 10^, and 
Viooo was taken as an approximation to the limit random variable V. We 
arrived at this choice using extensive exploratory analysis comparing Flood 
and Vn- The comparisons involved studying the empirical cumulative distri- 
bution functions and two-sample comparisons using Kolmogrov-Smirnov tests 
between Fiooo and other values of Vn- Specifically, based on 10^ samples and 
two-sample Kolmogrov-Smirnov tests, there were no statistically significant 
differences between Viooo and V200O) F500O) and Fio,ooo (with p- values > 0.185). 

Table 4.1 summarizes the probabilities of ruin (with M = 0) and the lower 
and upper bounds of the 95% confidence intervals (LCL, UCL) based on 10^ 
simulations. The confidence intervals in this example and other examples in 
this section are based on the simulations; that is, the lower 2.5% and upper 
97.5% quantiles of the simulated values of P{V > u). We also evaluated the 
true constant C of tail decay in Theorem 2.1 and the relative error (RE). Even 
in the extreme tail — far below the probabilities of practical interest in this 
problem — the algorithm works effectively and is clearly seen to have bounded 
relative error. For comparison, we also present here the crude Monte Carlo 
estimates of the probabilities of ruin based on 5 x 10^ realizations of V2ooo- 
We observe that for small values of u, the importance sampling estimate 
and the crude Monte Carlo estimates are close, which provides an empirical 
validation of the algorithm for small values of u. Admittedly, the value of the 
crude estimate for u = 10^ is questionable. Finally, since all of the conditions 
of Theorems 2.3 and 2.5 are satisfied, our algorithm is optimal and efficient, 
as can also be seen from the values of the relative error. 

Table 4.1. Importance sampling estimation for the ruin probability with 
investments obtained using M = 0. 



u 


P{V > u) 


LCL 


UCL 


C 


RE 


Crude Est. 


l.Oe+01 


5.86e-02 


5.65e-02 


6.07e-02 


2.33e-01 


1.84e+01 


5.73C-02 


2.0e+01 


3.66C-02 


3.52e-02 


3.81e-02 


2.21e-01 


1.98e+01 


3.54C-02 


l.Oe+02 


1.33C-02 


1.28C-02 


1.39e-02 


2.11e-01 


2.12C+01 


1.29C-02 


5.0e+02 


4.95e-03 


4.74e-03 


5.15e-03 


2.06e-01 


2.09C+01 


4.85e-03 


l.Oe+03 


3.27e-03 


3.14e-03 


3.41e-03 


2.07e-01 


2.12e+01 


3.21C-03 


5.0e+03 


1.25e-03 


1.19e-03 


1.30e-03 


2.06e-01 


2.18e+01 


1.23e-03 


l.Oe+04 


8.13C-04 


7.78C-04 


8.49e-04 


2.04e-01 


2.24e+01 


8.01e-04 


5.0e+04 


3.06C-04 


2.93e-04 


3.20e-04 


2.02e-01 


2.22e+01 


3.27C-04 


l.Oe+05 


1.98e-04 


1.90e-04 


2.07e-04 


1.98e-01 


2.16e+01 


2.10e-04 



4.3. Perpetuity sequences. A similar mathematical problem arises when 
one estimates the tail of a perpetuity sequence {Cn} as defined in (4.4), 
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namely 

(4.9) £* := Bo + AqBi + A0A1B2 + ■■■ . 

Note that by choosing Ai = R^ and Lj = Bi (and Bq = 0), we are back in 
the setting of the previous example (except that we now consider lim„_!.oo Cn 
rather than C := (sup^g^ £„) V 0, which we considered before). This leads 
to the SFPE 

(4.10) C=AC + B. 

If S > 0, then this is a special case of (4.5). 

The sequence (4.9) is of importance in life insurance mathematics. In that 
context, {An} denotes the discounted financial returns due, e.g., to inflation, 
while {Bn} denotes the future financial obligations of the company, generally 
taken to be positive. Then it is of interest to study P (£* > n) as u — )■ 00, 
i.e., the stationary tail on this sequence. 

Note that (4.9) can be viewed a backward sequence generated by the SFPE 
(4.10), while we simulate the forward sequence generated by the same SFPE. 
These two limiting distributions must be the same due to Letac's principle 
(c/. Letac (1986), CoUamore and Vidyashankar (2011), Lemma 2.1). 

4.4. The ARCH(l) process. We now estimate the tail probabilities of the 
ARCH(l) financial process, which were also studied in Blanchet et al. (2011). 
Originally introduced by Engle (1982), this process models the squared re- 
turns on an asset using the recurrence equation 

i?„ = (a + b-R„_i) Cn = AnRn^i + Bn, n = 1, 2, . . . , 

where An = ^Cn ^'^d Bn = aC,n ^^^ {Cn} is an i.i.d. Gaussian sequence. These 
assumptions imply E[logyl] < 00, and the additional assumptions of our 
theorems are easily seen to be satisfied. Setting Vn = -R^, we obtain that 

V := lim„_!.oo Vn satisfies the SFPE V = AV + B, which is the same as 
(4.10). Again, it is of interest to determine the tail behavior of {Vn} under 
stationarity, that is, P (V > u) for large values of u. Next we describe how 
our algorithm can be implemented to estimate these probabilities. 

4.4.1. Implementation. We set b = 4/5 and considered the following val- 
ues for a: 1.9 x 10^^, 1, and 2. It can be shown that 

(26)"r(« + l/2) 

^f^"J- — f(T7^^ — • 
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We solved the equation E [An] = 1 using Gauss-Hermite quadrature to obtain 
^ = 1.3438. Under the .^-shifted measure, An = bXn and Bn = aXn where 
Xn ~ r(^ + 1/2,2). Using the formulae in (4.8) for M, we obtained (upon 
taking the limit as /3 — t- and using the Taylor approximation r(/3 + 1/2) = 
r(l/2)+/3r'(l/2)+0(/32)) that Mq = 0.362, 0.724, 6.879x10-^ when a = 1, 2, 
and 1.9 X 10^^, respectively. Next we observed that M_g does not exist, 
since A(— ^) = oo in this example. However, using the numerical procedure 
at the beginning of this section, we deduced that, in fact, M_^ = 0. Thus 
we set M = Mq. We also computed Mq using the algorithm described at the 
beginning of this section and obtained similar estimates. 

Table 4.2 summarizes the simulation results for tail probabilities of the 
ARCH(l) process based on 10^ simulations. We notice a substantial agree- 
ment between the crude Monte Carlo estimates and those produced by our 
algorithm for small values of u. More importantly, we observe that the rela- 
tive error remains bounded in all of the cases considered, while the simulation 
results in the paper of Blanchet, Lam, and Zwart (2011) show that the relative 
error based on their algorithm increases as the parameter u — ;■ oo. 

4.5. The GARCH(1,1) process. A variant of the last example is the so- 
called GARCH(1,1) financial process introduced by Bollerslev (1986). Here 
the logarithmic returns on an asset are modeled as i?„ = anQn-, where o"„ 
denotes the stochastic volatility at time n and {Cn} is i.i.d. Gaussian. It is 
assumed that the squared volatility satisfies the recurrence equation 

(4.11) (T^ = Oo + &icr^_i + aii?^_i, n = l,2, 

Then V := lim„_!.oo o'n satisfies the SFPE 

(4.12) V = AV + B, where A = (6i + aiCi) and 5 = cq. 

To study the large exceedances of {i?„} under stationarity, it is of primary 
interest to determine P (U > u) as u — )■ oo. Let jX denote the distribution of 
A. Then it is easy to see that 



fi{{-oo,x\) 




for X < 6i, 
for X > bi, 



where ^{■) is the distribution function of a standard Gaussian random vari- 
able. Hence, under the .^-shifted measure, the random variable A has the 
probability law 

fi^{E) := / y^dii{y), E G B{R). 
Je 
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Table 4.2. Importance sampling estimation for the tail probability of 
ARCH(l) financial process with a = 1, 2, 1.9 x 10^^. 



u 


P{V > u) 


LCL 


UCL 


C 


RE 


Crude Est. 


l.Oe+01 


7.73C-02 


7.64e-02 


a== 1 
7.83e-02 


1.71C+00 


6.21C+00 


7.75C-02 


2.0e+01 


3.43e-02 


3.35e-02 


3.51e-02 


1.92e+00 


1.18C+01 


3.43e-02 


l.Oc+02 


4.34e-03 


4.23e-03 


4.45e-03 


2.11e+00 


1.29e+01 


4.28e-03 


5.0C+02 


5.07e-04 


4.96e-04 


5.18e-04 


2.15C+00 


1.13e+01 


5.21e-04 


l.Oe+03 


2.04e-04 


1.99e-04 


2.09e-04 


2.20e+00 


1.28C+01 


2.07e-04 


5.0e+03 


2.32e-05 


2.28e-05 


2.36e-05 


2.17e+00 


8.08e+00 


2.00e-05 


l.Oe+04 


9.00e-06 


8.88e-06 


9.12e-06 


2.14e+00 


6.83e+00 


9.00e-06 


5.0e+04 


1.07e-06 


1.05e-06 


l.lOe-06 


2.21e+00 


1.27C+01 


2.00e-06 


l.Oe+05 


4.11C-07 


4.04C-07 


4.18C-07 

a = 2 

1.64e-01 


2.15C+00 


8.51C+00 


NA 


l.Oe+01 


1.62C-01 


1.60e-01 


3.57e+00 


5.99e+00 


1.62C-01 


2.0e+01 


7.73C-02 


7.64e-02 


7.83e-02 


4.33e+00 


6.21e+00 


7.78e-02 


l.Oe+02 


1.08C-02 


1.05C-02 


l.llc-02 


5.25C+00 


1.34e+01 


1.06e-02 


5.0e+02 


1.28e-03 


1.25e-03 


1.31e-03 


5.43e+00 


1.14e+01 


1.33e-03 


l.Oe+03 


5.07e-04 


4.96e-04 


5.18e-04 


5.45e+00 


1.13e+01 


5.44e-04 


5.0e+03 


5.96e-05 


5.81e-05 


6.11e-05 


5.57C+00 


1.27e+01 


7.70e-05 


l.Oe+04 


2.32e-05 


2.28e-05 


2.36e-05 


5.51e+00 


8.08e+00 


3.50e-05 


5.0e+04 


2.64e-06 


2.60e-06 


2.68e-06 


5.44e+00 


7.60e+00 


3.00e-06 


l.Oe+05 


1.07e-06 


1.05e-06 


l.lOc-06 


5.61C+00 


1.27e+01 


l.OOe-06 






( 


a = 1.9 X 10-5 






l.Oe+01 


4.45e-08 


4.38e-08 


4.52C-08 


9.82e-07 


8.38e+00 


NA 


2.0e+01 


1.75e-08 


1.72e-08 


1.78e-08 


9.80e-07 


l.OOe+01 


NA 


l.Oe+02 


2.02e-09 


1.98e-09 


2.05C-09 


9.82C-07 


9.29e+00 


NA 


5.0e+02 


2.66e-10 


1.99e-10 


3.32e-10 


1.13C-06 


1.27e+02 


NA 


l.Oe+03 


9.59e-ll 


8.77e-ll 


1.04e-10 


1.03C-06 


4.38e+01 


NA 


5.0e+03 


1.04C-11 


1.02e-ll 


1.06e-ll 


9.75e-07 


l.Ole+01 


NA 


l.Oe+04 


4.15e-12 


4.05e-12 


4.26e-12 


9.85e-07 


1.32e+01 


NA 


5.0e+04 


4.78e-13 


4.66e-13 


4.91e-13 


9.86e-07 


1.34e+01 


NA 


l.Oe+05 


1.91C-13 


1.83C-13 


1.99C-13 


l.OOc-06 


2.19C+01 


NA 



Unlike in the ARCH(l) example, there is no closed form expression for fi^, so 
generating data from this distribution is more difficult. However, generating 
data from the shifted empirical distribution is easy, since it amounts to gen- 
erating data from a multinomial distribution. Let fij3^k denote the empirical 
distribution based on a sample of k observations from ^u^ for /3 = and P = ^] 
thus. 



A«,fc(-E^) := 



fRV^dfJ-oAv)' 



EeB{l 



In our simulations, we chose k = 10^ and a^ = 10 ^, oi = 0.11, and 6i 
0.88. Also, solving the equation E[A^] = 1 numerically, we obtained ^ 
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1.838. Furthermore, the minimum of ||A||q is 0.98. Hence Mq = SOao- Under 
the (— i^)-shifted measure, the minimum of ||^||a is and thus M_^ = ao- 
Therefore M = SOoo- The same value of M was obtained using the procedure 
described at the beginning of this section. 

Results for the estimation of the tail probabilities are summarized in Table 
4.3. As in the previous examples, the algorithm works effectively and is clearly 
seen to have bounded relative error. In addition, we considered the value 
ao = 0.01. In all of these simulation results, we see once again that the 
importance sampling estimate agrees very well with the crude Monte Carlo 
estimate for small values of u (where it is sensible to make this comparison) . 

Table 4.3. Importance sampling Estimation for the tail probability of 
GARCH(1,1) financial process. 



P{V > u) LCL 



UCL 



C 



R.E. 



l.Oc+01 


3.61e-12 


3.42e-12 


3.80e-12 


2.49e-10 


2.67C+01 


2.0C+01 


1.03e-12 


9.89e-13 


1.07e-12 


2.54e-10 


1.95C+01 


l.Oe+02 


5.24e-14 


5.05C-14 


5.43e-14 


2.49e-10 


1.88e+01 


5.0e+02 


2.58e-15 


2.52e-15 


2.64e-15 


2.36e-10 


1.17e+01 


l.Oc+03 


7.75e-16 


7.38C-16 


8.11e-16 


2.53e-10 


2.40C+01 


5.0e+03 


3.96e-17 


3.85e-17 


4.06C-17 


2.49C-10 


1.37e+01 


l.Oe+04 


1.09e-17 


1.07e-17 


1.12e-17 


2.46C-10 


1.28e+01 


5.0C+04 


5.78e-19 


5.59e-19 


5.97e-19 


2.51C-10 


1.67C+01 


l.Oc+05 


1.56e-19 


1.53e-19 


1.60C-19 


2.43C-10 


1.28C+01 



We end this section with a brief comparison of our work with that of 
Blanchet, Lam, and Zwart (2011). As shown above, our method works effec- 
tively and is comparatively easy to implement. The complex nature of their 
algorithm can be attributed to their alternative approach, which relies on us- 
ing complex asymptotic estimates. Furthermore, a critical issue with their al- 
gorithm is that it depends on several parameters, and effectively tuning them 
to produce optimal results seems to be challenging even when specialized to 
the simpler case of perpetuities without invoking Markov state-dependence. 
Additionally, their results on crude Monte Carlo — even for small values of u 
corresponding to A = 0.1 and 0.5 in their paper — are quite different from 
ours. We emphasize that our crude Monte Carlo estimates are based on five 
million simulations and on the approximation of the limit variable V by V2ooo- 
The latter approximation was based on extensive exploratory analysis. We 
are not able to explain this discrepancy with their results (as the implemen- 
tation details are not available in their manuscript). However, we emphasize 
that our crude Monte Carlo estimates and importance sampling estimates are 
quite close to one another for moderate values of u, as one should expect. Fi- 
nally, it is unclear from their numerical results whether their state-dependent 



IMPORTANCE SAMPLING FOR SFPE 35 

sampler yields unbiased point estimates or bounded relative error. 

5. Proofs of results concerning running time of the algorithm. 

We now turn to the proof of Theorem 2.4, which will rely on ideas from 
nonlinear renewal theory and Markov chain theory. In the proof, we will 
utilize the following: 

Lemma 5.1. There exist positive constants M and j3 and a constant p G 
(0, 1) such that 

(5.1) 'E.^[h{Vn)\Vn-l] < ph{Vn-i) on {K-1>M}, 

where h{x) := x~'^l|j.>i} + l{x<i}- 

Proof. Without loss of generality, assume that Vn-i = v > 1. Then by 
strong Markov property, 

Eg [h{Vn)\Vn-i = v] = B^[V^~H{v,>i}\Vo = v]+P^ {Vi < l|yo = v) . 

By assumption {H4), the second term on the right-hand side of the above 
expression is o(v''^). As for the first term, notice that it can be expressed 
as v^F,^ [{Ai uia.-K{v^^Di, 1) + v^^Bi)^^l^y^-^iy\Vo = v]. Using the bounded- 
ness of {(^1 m.a:xL{v~^Di, 1) + v^^Bi)~^lry^yij\Vo = w} in v, it follows that 
Ei^[{AiTaax{v^^Di,l) + v^^Bi)~/^l^Y^yij\Vo = v] converges as w — )■ 00 to 

Eg [^"^ ] = A(^ - /?) < 1 if < /3 < ^. Thus, choosing /3 = e G (0, 0, where e 
is given as in (-^4), we obtain that the lemma holds for any p = (Eg [A^"" ] , l) 
and M < 00 sufficiently large. □ 

Proof of Theorem 2.4 (i). Let M be given as in Lemma 5.1. Without 
loss of generality, we may assume that M > max(M, 1). Let L = sup |n G 
Z+ : Vn G (— oo, M]}. Then it follows directly from the definitions that K < L 
on {K < oo}. Thus it is sufficient to verify that Eg [L] < oo. 

To this end, introduce two sequences of random times, as follows. Begin by 
setting iTo = and /Cq = 0. Then for each i G Z+, set 

/Ci = inf{n > Ji^i :Vn> M} and Ji = inf{n > Id : Vn € (-oo, M]}. 

Intuitively, {/Ci} denotes the successive times that the process escapes from 
the interval (— oo, M], while {J^i} denotes the successive times that the process 
subsequently returns to (— oo, M], where < ICi < Ji < . . . . For any integer 
i, let Ki := K-i — /Cj_i denote the ith inter-escape time from the interval 
(— oo, M]. Finally, let 01 denote the total number of cycles that exit (— oo, M] 



i<E 
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and subsequently return to (— oo,M]. Then it follows from these definitions 
that 

Observe that E^ [9T] < oo. Indeed, when the process {Vn} escapes from 
(— oo,M], the probability that it ever returns to this interval is bounded 
above by this probability conditioned on the starting state Vq = M. But by 
the transience of {Vn} in Lemma 2.1 (ii), that probability is less than one. 
Thus 

sup P^ [Vn G (— oo,M], some n> l\Vo = v^ < p < 1, 

v>M 

and consequently E^ [9T] < Xl^iP" < °^- 

It remains to show that Eg [kj] < oo, uniformly in the starting state 
^Ki-i £ {M, oo]. But note that the random variable Ki can be divided into two 
parts; first, the sojourn time that the process {Vn} spends in (M, oo) prior to 
returning to ( — oo, M] and, second, the sojourn time in the interval (— oo, M] 
prior to exiting again. Now if K denotes the first return time to (— oo,M], 
then by Lemma 5.1, 



-^l{/^<oo} 



Vo 



< Q for some finite constant 0, uniformly in 



Hence Eg 

V > M. 

Thus, to establish the lemma, it is sufficient to show that Eg [iV Vq = f] < 
oo, uniformly in f G (— oo,M], where N denotes the total number of visits 
of {Vn} to (— oo, M]. To this end, observe that [—M,M] is petite (CoUamore 
and Vidyashankar (2011), Lemma 5.1). Moreover, it is easy to verify that 
(— oo, —M) is also petite for sufficiently large M. Indeed, under Letac's Model 
E, we have that m{^^t^_fj^P^{m.ax{Di,Vo) = Di\Vo = v) > 0, and in 
that case Vi = AiDi + Bi. Thus, the transition kernel of {V^} satisfies a 
minorization with small set (— oo, — M). Consequently (— oo, M] is petite and 
hence uniformly transient (Meyn and Tweedie (1993), Theorem 8.3.5, and the 
transience of {V^}). We conclude that Eg [iV] < oo, uniformly in the initial 
state Vo G (— oo,M], completing the proof. □ 

We now turn to the proof of Theorem 2.4 (ii). We begin with a slight 
variant of Lemma 4.1 in Collamore and Vidyashankar (2011). 

Lemma 5.2. Assume Letac's Model E, and suppose that {Hi), {H2), and 
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(Hs) are satisfied. Then 

(5.2) limpJ^yy Tu<K]=Pi:(v>y) 

for some random variable V . The distribution of this random variable V is 
independent of the initial distribution of Vq and is described as follows. If 
A is a typical ladder height of the process Sn = X^ILilogAj in the ^-shifted 
measure, then 

(5.3) P^ ( log y > y) = ^jj^ I P^ {A' >z)dz, for all y>0. 
Proof. Note by definition of Z„ that 

Vn = {Ao ■ ■ ■ An) {Znl{Z„yO} + Z„l{^^<o}) . 

Now it is shown in Lemma 5.2 of CoUamore and Vidyashankar (2011) that 
Vn'\ oo w.p. 1 as n — 7- oo under the measure ^g. Thus Tu < oo a.s., and at this 
exceedance time, we have that l|r^<o} = 0. Consequently, setting Xn = Zn 
on {Zn > 0} and Xn = 1 otherwise, we have that Vt^ = V^ , where 

K = (Ao---A„)X„, n = 0,l,.... 

Since Xn is positive for all n, we can introduce the perturbed random walk 

(5.4) \ogV;, = Sn + 5n, n = 0,l,..., 

where Sn := YJi=i fog-4i and 5n ■= logX„. 

Note that Sn has a positive drift under ^u^-measure and that {{log Ai,5i) : 
i = 1, . . . ,n} is independent of logylj for all j > n. Thus classical nonlinear 
renewal theory can be applied to the sequence {logV^}. To do so, we need to 
verify that the sequence {Sn} is slowly changing. But by Lemma 3.1, we have 
that lf^z„<o} ~^ a.s. and that (5„ := logXn converges to a proper random 
variable. Hence (5„/n — )■ a.s., and consequently {5n} is slowly changing. 
Hence it follows from Theorem 4.2 of Woodroofe (1982) that 

(5.5) -^ ^V as u^oo. 

u 

Note that this result holds independent of the initial distribution of Vq G C. 
Next observe that for any y > 1, 

(5.6) V^{^>y^=V^{^>y-T^< K^ +P^ (^Ik > y ; T„ > K 
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By the strong Markov property, 



(5.7) p^ ( ^ > y 

u 



T^>K]=Pi:{^^>y 



V^ 



u 



Vq ^ u 



where ly is the distribution of V^ conditional on {T„ > K}. Note by definition 
of K that y^ G C. Hence it follows by (5.5) and (5.7) that 

/V' 
(5.8) lim P^ -^ > y 



-^i — ^yTu>K\=P^[V>y], y > I. 
Consequently by (5.5) and (5.6), 

(5.9) JUm P^ {^ > y|T„ < i^") = P^ (y > y) , y > 1, 

provided that lim inf ^^oo Pg (Tu < K) > 0. Finally, to verify that this last 
condition is fulfilled, observe that P^ (T„ < K) increases to P^ (K = oo) as 
u — )■ oo. Moreover, by Lemma 5.2 of Collamore and Vidyashankar (2011), we 
have that T^ t co w.p. 1 as n — )■ oo in the ^-shifted measure. Thus, starting 
from the stationary measure 7 on C, P^ {K = 00) is strictly positive over any 
given cycle. Thus P^ (T^ < K) — )• as n — ?■ 00, for some positive constant 0. 

D 

Proof of Theorem 2.4 (ii). Let V^ be defined as in the proof of Lemma 
5.2, and observe from the proof of this lemma that logl/^' = Sn + Sn, where 
{6n} is slowly changing and Sn = Y17=i ^^gAi. 

Set T^ = inf{n : V^ > u}. Then it follows by Lemma 9.13 of Siegmund 
(1985) that 

T' 1 

(5.10) -—^ > .,,^, in probability 

with respect to the measure /i^, since A'(,^) = E^[log^]. Now if T„ < K, then 
the process {Vn} never returns to the set (—00, 0], and hence we have in that 
case that T^ = T^. Moreover, since Vnt "^^ w.p. 1 as n — )■ 00 in the ^-shifted 
measure, P^ (r„ < K) — )• O as w — )■ 00, for some positive constant 0. Thus 
it follows from (5.10) that, conditional on {T„ < K}, (Tu/logu) -^ (A'(C))"-^ 
in probability. 

To show that convergence in probability implies convergence in expec- 
tation, it suffices to show that the sequence {Tu/logu} is uniformly inte- 
grable. Let M be given as in Lemma 5.1, and first suppose that M < M and 
supp (Vn) C [— M, 00) for all n. Then in this case, conditional on {T^ < K}, 

Tu> n =^ Vi G (M, ti), i = 1, . . . , n. 
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Now apply Lemma 5.1. Iterating (5.1) — as in the proof of Lemma 3.3 — we 
obtain that (3.29) holds with p in place oi pp. Consequently, using the explicit 
form of the function h in Lemma 5.1, we obtain that 

(5.11) VdTu> n\ Tu<K)< L^^^^——\ p-u^, for ah n, 

where /3 is given as in Lemma 5.1. Now P^ (T„ < K) J, > as u — )■ oo. 

Hence, letting El" [•] denote the expectation conditional on {T„ < i^}, we 
conclude that 



(5.12) eJ"^ 



T T 

\ !i - ^ 

log u log u 



for some finite constant 0, and for sufficiently large 77, the right-hand side 
converges to zero as u — )■ 00. Hence {r„/logn} is uniformly integrable. 

If the assumptions at the beginning of the previous paragraph are not 
satisfied, then write T„ = L + (T„ — L), where L is the last exit time from the 
interval (—00, M], as defined in the proof of Theorem 2.4 (i). Then (T„ — L) 
describes the length of the last excursion to level u after exiting (— oo,M] 
forever. By a repetition of the argument just given, we obtain that (5.11) 
holds with {Tu — L) in place of T„, and hence (5.12) also holds with (T^ — L) 
in place of T^. This implies {(T„ — L)/\ogu} is uniformly integrable. Next 
observe by the proof of Theorem 2.4 (i) that E^ [L/logu] J, as u — t- 00. The 
result follows. □ 

In Theorem 2.4 (ii), we showed that 

(5.13) E^[r„|T„<K]-^ as n ^ 00, 

which can be expected due to the fact that A'(i^) = E^ [log A]. Specifically, if 
we were to replace {Vn} with the multiplicative process {Ai ■ ■ ■ An} and define 
Tu = inf{n : {Ai---An) > u}, then by an elementary renewal argument 
applied to Sn := ^l=i^ogAi, we would obtain that (5.13) holds with the 
left-hand side replaced with E [T^] . 

Next we turn to the proof of Theorem 2.4 (iii), which provides a similar 
asymptotic estimate for the return time to the set C, assuming that T^ < K 
and starting our analysis from the state Vt„ > u. Now if T^ < K, then 
starting at time T„, we have Vr^ ~ u and we would like to study its return 
time to C = [—M,M]. Again, if we replace {Vn} with the multiplicative 
process {{Ai • • • An)u} and define K = inf{n : {Ai • • • An)u < M}, then an 
elementary renewal argument yields 

_, r ~ -, log u 
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where we have used that E [log A] = A'(0) < 0. This heuristic argument serves 
as motivation for Theorem 2.4 (iii). 

We will first establish Theorem 2.4 (iii) conditional on the event that 
(Vb/u) = V > 1, and later remove this assumption. Thus, we assume for 
the moment that the process starts at a level vu where f > 1, so the dual 
measure for the process agrees with its initial measure. For this process, define 

L(z) = inf{n : \Vn\ < z}, for any z >0. 

Lemma 5.3. Let (Vq/u) = v > 1 andt G (0, 1). Then under the conditions 
of Theorem 2.4, 



(5.14) lim E 



u—>-oo 



lU) 



Vo 

V 



logu I ^ ^ u J |A'(0)|' 



1 



Proof. For notational simplicity, we will suppress the conditioning on 
(Vq/u) = f in the proof of the lemma. We begin by establishing an upper 
bound. First observe that, as in (3.27), 

(5.15) jK\^^^^{A^m + \B^\) 



Hence 

(5.16) log(j^^^ < log (^An + u-'{An\Dn\ + \Bn\)), n<L{u'). 

This shows that we may bound the process {\Vn\ '■ n < L{u^)} by a classical 
random walk. More precisely, define 

n 

5i") := J2 ^t^ ^here xf") := log (^A^ + u'' (^.| A| + m) ) . 

Then iterating (5.16) and using that [Vq/u) = v, we obtain 

(5.17) log IKI - log{vu) < 5i"\ for all n < L{u^). 

Now let L„(u*) = inf {n : S^^ < -(1 - t) logu- logt;}. Then L(n*) < Lu{u^) 
for all u. 

Since {Sn j is a classical random walk and Lu{u^) is a stopping time, it 
follows by Wald's identity that ^[Si^^^f^] = E[x|"^]E[L„(u*)] . Thus letting 

On ■■= l-^Z^K) ~ (1 ~ *) logM - logvl 
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denote the overjump of the random walk {S*™ } over its boundary (which in 
this case is taken to be the level (1 — t) logu + logw), we obtain 



(5.18) L(u*) < 



(1 - t) logu + log t; + E[0„] 



|E[X1 



Hi 



Since E[X}" ] — )• A'(0) as n — )■ oo, the required upper bound will be estab- 
lished once we show that 

(5.19) lim E[OJ =0. 

w^oo log U 

To establish (5.19) note, as in the proof of Lorden's inequality (Asmussen 
(2003), Proposition V.6.1), that E [0„] < E [Y^] /E[y„], where Yu has the 

negative ladder height distribution of the process {Sn }■ Next observe by 
Corollary VIII.4.4 of Asmussen (2003) that 

(5.20) E [Yu] = m(,i)e®" -^B[Y] as u ^ oo, 

where Y has the negative ladder height distribution of the process {Sn}, 

oo ^ 

m(f):=|E[x(«)]|,j = l,2,..., and 6„ := ^ -P (^i") > o) . 

n=l 

We observe that ©„ is the so-called Spitzer series. Similarly, an easy calcula- 
tion (c/. Siegmund (1985), p. 176) yields 



oo ^ p , - 

(5.21) E [y„^] = mf e^" -2m«e^" J] ^E (si^A ^ E [Y 



oo 

n — )■ OO. 



n 

n=l 



Since E [(log Af] < oo ^ E [Y^] < oo for j = 1, 2, it follows that E [On] -^ 
E [y^] /E [Y], i.e. to a finite constant, which implies (5.19). Thus (5.14) holds 
as an upper bound. 

To establish a corresponding lower bound, fix s G (t, 1) and define 

L{u^) = inf {n : Sn < —{I — s) logu — logv}. 

(Essentially, L{z) is defined in the same way as Lu{z), except that we now 
substitute the sequence {Sn} in place of {Sn }.) Next observe that Vn > 
^nKi-i — \Bn\ provided that Vn-i > 0, and iterating this equation yields 

n n 

(5.22) Vn>(Ai--- An) Vo - W, where W := lim V TT AABA. 

?i— >-oo ^ — ' J--'- 
1=1 j=i+l 
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Now from the very definition of L and the assumption that {Vq/u) = v, we 
obtain that 

L{u'') > n ^=> (Ai • • • Ak)Vo > u", for aU k < n. 

But by (5.22), (yli • • • Ak)Vo > u' ^ Vk > u^ on {W < {u' - u*)}. Thus for 
all n, L{u^) > n =^ L{u^) > n on {W < {u^ — n*)}, and consequently 

(5.23) E [L(u*)] > E \l{u'); W < {u' - u^) 

To study the expectation on the right-hand side, recall that W satisfies 

the SFPE W = AW + \B\, and hence by Theorem 2.1 of Collamore and 
Vidyashankar (2011), 

(5.24) ¥[W >u' - n*) ~ Cu-'^ as n ^ oo. 

Next observe that L{u^) is, by definition, the time required for the negative- 
drift random walk {S'n+logf } to reach the level — (1 — s) logu. But by Heyde's 
(1966) a.s. convergence theorem for renewal processes, 

L(u^) (1-s) 

(5.25) ^1 > , . ,, , , a.s. as u — 7- oo, 

^ ' logu |A'(0)| 

where in the denominator of this last expression, we have observed that 
E [log^] = A'(0) < 0. It foUows from (5.25) that for any e > 0, 

( Liu*) ^ , ^^ , 1-s 

(5.26) lim P -^ — ^ ^ (r - e,r + e) = 0, where r :- 



logu ^ ' / l^'(0)l 

Substituting (5.24) and (5.26) into (5.24) and letting e — )■ 0, we obtain that 

(5.27) liminf r^E \L (u*)1 > ,^~tn 
^ ' M^oo logu L V ;j - |A'(0)|' 

and the required lower bound is obtained by letting s \^t. D 

Lemma 5.4. Under the conditions of Theorem 2.4, 

(5.28) lim i lim sup E \l(M) - L (u*) 11=0. 

40 i u^oo logu L ' ' J J 
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Proof. Applying Lemma 3.3 with /? = 0, we obtain that the drift condi- 
tion (T>) holds, and hence for some a > 0, 

E[|V;|" I T4_i = w] < p|ti;|", for aU w^C, 

where p £ (0, 1), and C = [—M,M] (where M > Mq and Mq is given as in 
Lemma 3.3). Iterating this equation yields 



(5.29) 



E[l{L(M)>n}\Vo = w] <P"(^) 



for all n. 



We now apply this result to obtain an estimate for L(M) — L{u*'). Recall 
that L{M) — L{u^) measures the length of time required for the process, 
beginning at level V£,(„t), to enter the set C = [— M, M]. Thus by definition, 
V/;,(ut) < u^ ■ Hence we are in the setting of (5.29) with w < u*, and with 
L{M) replaced by L{M) — L{v^). Using the strong Markov property and 
(5.29), we then obtain that 



(5.30) 



P L(M)-L(n*) >n <p 



M 



t \ a 



for all n. 



Set Jt{u) = L{M) — L{u*) and t' = to/(— logp). It follows from the previous 
equation (upon summing over all n > t' log u) that 



(5.31) 

Then 

(5.32) 



E 



<^t (")l{Jt («)>*' log «} 



< 



pt'logu f yt\ a 



1-/9 V^-^ 



1 



(1 - p)M^ 



lim sup 



logTi 



-E 



L{M) - L (u*) 



<t'. 



Since i' J, as t J, 0, we conclude (5.28). 



D 



Proof of Theorem 2.4 (iii). Set 

1 



HJv) 



logii 



-E 



L{M) 



Vq 



u 



Then it follows from Lemmas 5.3 and 5.4 that 

lim Huiv) = — 



(5.33) 



A'(0)|' 



Let flu, A denote the probability laws of the random variables Vr^/u, V, 
respectively, as given in Lemma 5.2. Then using the strong Markov property. 
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it follows that L(M), conditional on Vq/u ^ fiu, is equal in distribution to 
K — Tu, conditional on {T„ < K}. Thus it is sufficient to verify that 



(5.34) lim 



f 



logu 



-E 



L{M) 



Vo 



u 



fJ-u 



lim / Hu{v)d(iu{v) 



1 



|A'(0)| 



This result will follow from (5.33), provided that we can show that the limit 
can be taken inside the integral in the above equation. 
To do so, express the inner quantity in (5.34) as 



(5.35) 



Hu{v)d {flu - A) iy) + / Hu{v)dfi,{v). 

v>0 Jv>0 



To deal with the first term, we begin by obtaining an upper bound for Hu{v). 
First note by a slight modification of (5.31) (with t = 1, 9 = —logp, Ji{u) 
replaced with L{M), and t' replaced with an arbitrary constant r) that 



(5.36) 



E 



LiM)l{^M)>r\ogu} |Vb| < u 



< 



-9 / 
I U 

l-~p\M 



for all r > and some a > 0. Now choose r > a/9. Then by the previous 
inequality, 



(5.37) 



E 



^Wl{L(M)>rlog«} I Vol <U 



<ei, 



for some finite constant ©i which is independent of u. Consequently 



(5.38) 



1 



\ogu 



-E 



L{M) |Vb| < u 



< r + 



logu 



Next, we extend this estimate to the case where the starting state satisfies 
{Vq/u) = w > 1. To this end, we write E [L(M)|(Vb/^) = i'] as a sum of two 
terms; first, the expected time for the process {Vn} to reach C = [—M.,M] 
starting from an initial state in [— n, w]; and second, the expected time to 
reach [— u,u] starting from the initial state (Vq/u) = v. Using the strong 
Markov property, this yields 



(5.39) E 



L{M) 



u 



< sup E 

«ie(Af,M] 



L{M) 



IVn 



w 



+ E 



Liu) 



u 



For the second term, we bound L{u) by the sojourn time for a classical random 
walk, as follows. Begin by observing that 



|K-i| > u 



\Vr,\ < IK 



n— 1 



Af, + 



Bn 
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Now taking logarithms, we see that E [L{u) {(Vq/u) = v] is bounded above by 
the length of time needed for the random walk 



'5i"^:=^r-i + logU„ + 



B„ 



u 



n 



1,2,. 



starting from Sq = log(fM), to reach the level logu. Denote this sojourn 
time by L*{u). Applying Lorden's inequality (Asmussen (2003), Proposition 
V.6.1) to {5n }, we obtain 

E[L*(n)] <Q2{u)\ogv + e^{u) 

for constants 02 ("u) — )■ m^ and 03 (n) — )■ m2/'mf, where mi and m2 are the 
first and second moments of the ladder height distribution for the sequence 
{log^i}; cf. the discussion following (5.19) above. (The required moment 
condition needed for this convergence to hold are satisfied due to the dif- 
ferentiability assumption in (Hi).) Finally, substituting this last bound and 
(5.38) into (5.39), we obtain that for some constant 0, uniformly in u > uq, 



(5.40) 



Hu{v) :-- 



1 



logu 



E 



L{M) 



Vo 



<0 + 



logw 
mi 



Returning to (5.35), using the above upper bound, we now show that 



(5.41) 



+ i5^)d(/i„-/i)(t;) 

v>0 \ "T-l 



as u 



oo. 



Since /i„ =^ jl, by Lemma 5.2, it is sufficient to show that J^^^logv djlu{v) is 
uniformly bounded in u, which would follow from the uniform integrability of 
{| log Vt„ — logu|}. To this end, we apply the Corollary to Theorem 2 of Lai 
and Siegmund (1979). First, write the process {Vn} as a nonlinear renewal 
process as in (5.4). Then apply Lemma 3.1 to obtain (in the same notation 
as in (5.4)) that 



(5.42) 



Sn < log V 



Bi 



1=0 



A, 



^{K>i} 



n < K. 



(This is the first equation in (3.1).) Next apply Jensen's inequality and Lemma 
3.1 (ii) to obtain that 



(5.43) 



e 



E, 



\^Tu - ^Tu-l\ l{T„<i^} 



< E. 



\og {Z(^^Y 



< CX). 
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Notice also that conditions (6)-(8) of Lai and Siegmund (1979) are satisfied 
(with a = 1). By (Hi), we have that E^ [(logA)^] < oo. By a sUght modifi- 
cation of (5.43) followed by an application of Chebyshev's inequality, Eq. (7) 
of their article follows. Moreover, we observe that Theorem 2 of their article 
is actually valid if their Eq. (8) is replaced by uniform continuity in proba- 
bility of {5n}, as given in Eq. (4.2) of Woodroofe (1982). The latter property 
holds, since (5„ converges w.p. 1 to a proper random variable. We conclude 
that all the conditions in Lai and Siegmund (1979) are satisfied, and hence 
{| log Vt„ — logu|} is uniformly integrable. Consequently (5.41) follows from 
the weak convergence /i„ =^ /i. 

Finally, applying the dominated convergence theorem to the second term 
in (5.35) and using (5.33), we conclude that (logu)~-'^E [L{M) \(yo/u) ~ /i] — ;■ 
1/|A'(0)|, as required. □ 

6. Proof of optimality. Next we turn to the proof of the optimality the- 
orem. We remark that a similar result was obtained in a different setting in 
CoUamore (2002), although the current proof incorporates new aspects, most 
notably, the possibility that the alternative algorithm be state-dependent. 

Proof of Theorem 2.4. We divide the proof into two parts, "asymptotic 
equivalence" and "uniqueness." 

Part I: Asymptotic equivalence. We begin by showing that if \^u{E;w,q) : 
E G i3(M^), tf G M, g G {0, 1}} is any family of measures within the class 93T, 
where 9JT is given as in (2.15), then 



3.1) 



lim inf 

u— >oo log U 



E,. 



(4^))^ 



> 



lim 

«^oo log U 



^£1] 



-2^; 



in other words, simulation under the dual measure used in our main algorithm 
is either asymptotically equivalent or preferable to any alternative in 9Jt. 
To establish (6.1), set 



H^{E;w,q) 



l^dE), 


E £B{R^), w£R, and q = 0; 


KE), 


E el3{R^), weR, and g = 1. 



(In this definition, recall that intuitively, w corresponds to the level of the 
process {logl4i_i/logn}, while q = 1 corresponds to the event that {Vn} has 
exceeded level u by the previous time.) 

First assume that i^ <C ^2)- Then by applying Lemma 4.2 twice, first to 
express the given expectation in terms of the original measure and then in 
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E,. 



(4'^))' 



E. 



Ej) 



d^ 



K\t^<k}\{[-£{Y,;W,,Q, 

K y 



where Yi := (log Ai,Bi,Di) is the driving sequence of the SFPE, Wi := 
Vi-i/u, and Qi is the indicator function on {Vj > u,for some j < i}. Using 
the Radon-Nykodym theorem to simphfy the last quantity in the previous 
display, we deduce that 



E„ 



[sir^f 



E^ 



D 



K 



^«1{T„<A-} 



dfi 






dv 



{Yi;Wi,Q^ 



Observe that by the definition of the dual measure, j^ = j^ when Qi = 
0, because in this case, the process {Vn} has not attained the level u by 
time i, and consequently the dual measure agrees with the ^-shifted measure. 
Similarly, ^rf^ = 1 when Qi = 1, because in that case, the process {Vn} has 
attained the level u by time i, in which case the dual measure agrees with the 
original measure. Thus, from the previous equation we obtain that 



(6.2) eJ(£:H)^ 



=E 



D 



dH , , \ T-r d\xi 






Now set 



Ui = log 



dv 
d^^ 



{Yi■,W^,Q^)] and i?„ = ^t/,. 
^ 1=1 



Substituting these definitions into (6.2) and noting that 3^(^) = e «^°s^i 
we obtain that 

(6.3) eJ(4'^))'1=E2, 



Nu^{T^<K} e 



-2^St^-Rk 



By Jensen's inequality and the observation that N^ > 1 on {T„ < K}, it 
follows that 



(6.4) 



E,, 



£. 



(uU2 



> pu exp {Ej, [-2^St^ - Rk\ T„ < K] }, 
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where 

Pu := P^ {Tu < K) — )■ as u — )■ oo. 

where Q := P^ {K = oo). Note that is positive due to the transience of 
{Vn} in its ^-shifted measure, as stated in Lemma 2.1 (ii). Consequently (6.4) 
yields 



1 



logii 



'H^2' 



1 



.5) liminf logE^ {8\^>y > - limsup Eg [2^St^ |T„ < K] 



M— >00 



logn 



limsup- E^[Rk\Tu<K] 

logu 



u—>-oo 



To identify the first term on the right-hand side of (6.5), note by Wald's 
identity that 

Eg [log^] Eg [T„ A K] = Eg [StA{t^<k}] + Eg [Sk1{t^>k}] , 

which yields 

Eg [Tu A K] Eg [Sk1{t^>k}] 



(6.6) Eg[5rjT„<K]=Eg[log^] 



Pg(T„<K) Pg(T„<K) • 



Now consider the first term on the right-hand side. Since Eg [T„ A K] 
Eg [Tu1{t^<^k}] + Eg [K1{t„>k}] , we obtain that 

Eg[T„AK] _ Eg [Kl{T„>x}] 

^'•'^ Pg (T„ < K) - ^^ [^" "^« < ^] + Pg(r„<i^) • 

By Theorem 2.4 (ii), 

(6.8) lim — ^Eg[r„|r„<K] 



logu^^ "' " ^ ^ A'(0' 
Moreover, as u — )■ oo, 

(6.9) Eg [i^l{T„>i^}] ^ Eg [Kl{i^<oo}] < oo, 

where finiteness on the right-hand side is obtained by Theorem 2.4 (i). Since 
Pg {Tu < i^) ^ Pg (-fC = oo) = e > 0, by Lemma 2.1 (ii), it follows from 
(6.7) and the previous two equations that 

u-^^logu\Pf{Tu<K)J A'(0 
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Next consider the second term on the right-hand side of (6.6). By Corohary 
3.1, 



E, 



e-«^^ 1{T^>K}] = E [1{T„>X}] ^ P (i^ < oo) as 



U —7- oo. 



Moreover, by the recurrence of the Markov chain {Vn} in its original measure, 

P (i^ < oo) = 1. Thus 



1 = hm Et 



e-«^-l{T„>/.} 



> hmsupE^ 



e«'^^'l{5K<o}l{T„>/^} 



Applying the inequality e^^ > 1 + ^x, a; > 0, to the expectation appearing on 
the right-hand side of the previous equation, we obtain 



3.11) 



lini sup 

«^oo log U 



E^ [Sk1{t^>k}] ) < 0. 



Finally, substituting (6.10) and (6.11) into (6.6), using that P^ (T„ < K) 
> 0, and recalling that E^ [log^] = A'(^), we conclude 



3.12) 



limsup E^ [5t„ \Tu <K]<1. 

u^oo log U 



We now turn to the second limit on the right-hand side of (6.5). Our 
objective is to show that 



(6.13) 



limsup Ej, [Rk \Tu < i^l < 0. 

u-¥oo log U 



To establish (6.13), first introduce the assumption that 

(A) log ( ) is bounded from below by a finite constant. 

This assumption will later be removed. 

Claim: Under (A), M.n '■= Rn^'S [1{t„<x}| 5^n] is a supermartingale. 
Proof of the Claim: Recall 

Rn = Yl Ui^ ^here U^ := log i^{Y^■, W,, Q,)j • 

Thus 



E® [7W„| 5^n-l] — -A^n-l+Ej) [C/nEx) [ l|7;^<i^} | S"™] | ^Jn-l] . 
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In the last term on the right-hand side, the inner conditional expectation 
is bounded above by one, and hence this entire term is bounded above by 
Ej) [Un I (Wn, Qn) ] , where we have used the strong Markov property to replace 
^n-i with {Vn-i-, l{T„<„-i}), or equivalently iWn-, Qn)- Now condition on the 
event that iWniQn) = {w,q) and apply Jensen's inequality to deduce that 

^^[Un\{Wn,Qn) = {w,q)]= / log ( -^ — {y]w,q)]dfi^{y]w,q) 



< log / diy{y;w,q) = 0, 

where, for notational simplicity, we have suppressed the dependence on (w, q) 
in the above integrals. Since this inequality holds for any (w, q), we conclude 
that 

(6.14) B^[Mn\dn-l]<Mn-l. 

Next, after using the inequality logx < x, observe that 

^^[U+\iWn,Qn) = iw,q)]< I dv{y-w,q)<l, 

JG+(v,q) 

where G+{v,q) := {y : log {-^{y;w,q)) > O} and [/+ = max(;7„,0). 

Moreover, since we are assuming that log ( j^ I is bounded from below by 
a constant, it follows that Un is also bounded from below by a constant, 
uniformly in n (where this constant is independent of ^n-i)- Thus, setting 
U^ = — min([/„, 0), we obtain that 

B^[U-\{Wn,Qn) = {w,q)] <@o, 

where Go is a finite constant independent of {w, q). Thus 



3.15) E \Mn-M 



hi-l 



d. 



n~l 



<eo + i. 



It follows that {Mn}nez+ lies in L^. By (6.14), we conclude that {Mn}nez+ 
is a supermartingale. □ 

Returning to the proof of the theorem, observe that since {Mn} is a su- 
permartingale, it follows from the optional sampling theorem that 

(6.16) Es [Rk1{t^<k}] = Ed [Mr] < 0, 
implying that 

(6.17) limsup Ex, [Rk1{t^<k}] < 0. 

M— 5>oo iOg U 
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Substituting (6.12) and (6.17) into (6.5), we conclude that 



(6.18) 



lim inf 



1 



-E,. 



(£,^^^) 



>-2e. 



w^oo log U 

This is the required result, since it is known by Theorem 2.3 that 

u^oo iOg U 

Thus we see that simulation under the dual measure is either asymptotically 
equivalent or optimal when compared with any other member in 931. 

If log ( j^) is not bounded from below by a constant, then in (6.2) we can 

replace u with a larger measure, namely v^*^' := v + e//^), where e > 0. Then 
(6.2) yields 



3.19) E^ 



{&^) 



>Ei, 



«:^w.<.,n(|J(>-.))nlS«:'^o.*: 



i=i 



Since log ( -^ — I > loge, the entire proof can be repeated without change, 
except that when applying Jensen's inequality we now obtain 

Es[C/„|(W„,Q„) = {w,q)] < log / diy^{y;w,q) = log(l + e) < e, 

and hence in this case the optional sampling theorem yields 
(6.20) Ej, [Rk1{t^<k}] < eEs, [K] . 

This last expectation can be decomposed into three terms, 

E^ [Kli^K<Tu}] + ^i [^«1{T„<X}] + Ed [{K - Tu) l{T„<i^}] • 
Hence, by applying Theorem 2.4 (i)-(iii), we obtain 



(6.21) 



\"^^^P bk^^ [^^1|^.<^}] ^ ^ ( A^ + ^ 



But e > was arbitrary. Thus we conclude that (6.18) holds in general. 

Part II: Uniqueness. It remains to show that strict inequality holds in 
(6.18) when the family {i^(-)} differs from the dual measure /Uj). First observe 
by Jensen's inequality that 

diy 



B^[Un\iWn,Qn) = iw,q)] 



log 



dfi's 
< log / du{y;w,q) = 



{y;w,q) dn^{y;w,q) 
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where, as before, we have suppressed the dependence on (w, q) in the above 
integrals. Moreover, for any given {w,q), equaUty holds in this inequality if 



and only if -j^ is equal to a constant a.s. Thus, ii u j^ /ij), then there exists 
a point [w, q) where 



(6.22) 



E® [Un\{Wn,Qn) = iw,q)] = -A, for some A > 0. 



Our next objective is to use the continuity condition (Co) to extend (6.22) 
from a point w G [0, 1] to an interval. The continuity assumption (Cq) states 
that 

Es, log (^{Yn, Wn, Qn)] \Wn = W 

is piecewise continuous as a function of w, for any fixed Qn G {0, 1}. We claim 
that, as a consequence. 



(6.23) 



E© 



U„ 



Wn = W 



E, 



log 



du 
dux, 



\^m l^nj Wn, 



Wn=W 



is also piecewise continuous as a function of w. To establish (6.23), recall that 
we are assuming ^ <^ u and u <C /Uji. (The latter condition will later be 
dropped.) Then by the Radon-Nikodym theorem. 



log 



dv) 



log 



dv 



dfi: 



D 



+ log 



dfi-s 



Since the expectation of the second quantity on the right-hand side is clearly 
piecewise continuous as a function of w, (6.23) follows immediately from the 
piecewise continuity in (Co). Using the piecewise continuity in (6.23), observe 
that (6.22) implies that this equation holds for all (w, q) & Gx {r}, where C is 
a neighborhood containing the original point w in (6.22) and where r G {0, 1}. 
We now show that by sharpening the estimate in Jensen's inequality on the 
set G X {r}, we obtain a strict inequality in (6.18). First assume that log (37^) 
is bounded from below by a constant, i.e., that condition (A) holds. Then 
by repeating our previous arguments, but using the sharper estimate (6.22) 
(when (w, q) G G x {r}) together with Jensen's inequality (when {w, q) ^ 
G X {r}), we obtain 

X; :=([/* + ... + [/*)Ex,[l|r„<^}|5„] 

is a supermartingale, where U* := Ui + Al^^/^^Qylj^Q^^j.y. Applying the 
optional sampling theorem to this process {7W*}, we deduce that 

(6.24) Ex, [Rk1{t^<k}] < -A {l{.=o}Ej, [of] + l|,=i}Es, [ol^] } , 
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where 

of := j;i|H.„eG} and 0« := ^ l{^„eG}; 

n=0 n=Tu+l 

that is, Ou denotes the occupation time the scaled process {logy„/logn} 
spends in the interval G during a trajectory starting at time and ending at 
time Tu-, while Ou denotes the occupation time that {\ogVn/\ogu} spends 
in the interval G during a trajectory starting at r„ and ending at time K. 

First assume that r = 0, so we are in the first scenario, where we would 
like to study the occupation measure Ou of the process {log Vn/ log u] in the 
set G prior to time Tu or, equivalently, the occupation measure of the process 
{Vn} in the set G{u) := [tt*, «*] prior to time Tu- First introduce a subinterval 
G'{u) := [u^ ,u*] strictly contained in [u^,u*]. Since {V^} is transient in its 
^-shifted measure, by Lemma 2.1 (ii), P^ (T^ < K) — t- P^ (K = oo) = > 
as n — )■ oo. Thus {Vn} enters the region [u^ ,oo) with a positive probability 
which tends to as it — )• oo. We claim that, in fact, {Vn} enters the smaller 
region G'{u) with a positive probability tending to as u — )■ oo. To establish 
this claim, observe by Lemma 5.2 that conditional on {T„ < K}, 

log Vti, =^ log V + log V as w — )■ oo 

for a proper random variable V supported on (0,oo). Consequently, setting 
viu] = n* , we obtain that 



3.25) P^ (^Vt„(„) G G'{u) and r,(„) < K] 

—?■ P^ (logy < {t' - s') logu and T^(u) < K) —^ @ as u —^ oo. 



Now consider the process {Vn} after it enters the region G'{u). In particular, 
suppose that this process starts with an initial state Vq G G'{u). Then our 
objective is to show that the expected number of visits of {Vn} to G{u) prior 
to termination at time Tu behaves, roughly speaking, like a constant multiple 
of logu as u — )■ oo. 

Fix 5 < min(s' — s,t — t'). Then observe by a law of large numbers argument 
that given 6, there exists a positive constant r] = i]{6) such that 

(6.26) Pg (I'S'nl > (5 logu, some n < rjlogu) — ;■ as u — )■ oo, 

where Sn = J2i=i^'^S^i- Next observe that if Vn-i > 0, then the equation 
Vn = An max (D„, Ki-i) + Bn yields the inequalities 

A-nVn-l — Bn < Vn < AnVn-l + -Bn, 
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where Bn := (yl„|-Dn| + |-Bn|). Iterating these equations yields 

V °° B 

(6.27) Vo-W < J- ^^— - <Vo + W, where W :=Y^ 



{Ai---An) - >^ ' ' ■ Z-.AI---A. 



i=l 



Hence, if Vq G G'{u) := [n'' , n* ] and y„ ^ G(u) := [u^,m*] for some n < 
rjlogu, then we must either have 

(6.28) {Ai---An)^[u'\u^], for some n < 7? log ti, 
or 

(6.29) W > min (u'' - u'+\ u^-^ - u^' 

But in the ^-shifted measure, M^ is a perpetuity sequence which converges 
to a proper random variable (c/. the proof of Lemma 5.2 in CoUamore and 
Vidyashankar (2011)). Hence for sufficiently small S, the probability in (6.29) 
tends to zero as m — )■ oo. Similarly, the probability in (6.28) also tends to zero 
as M — )• oo, by (6.26). Thus we conclude 

(6.30) P^{Vn ^ G{u), some n<r] log u\Vo £ G'{u)) ^ as u -> oo, 

uniformly in the initial state Vq G G'{u). Combining this result with (6.25) 
and invoking the strong Markov property, we obtain that 

(6.31) liminf ^WA > ri > 0. 

^ ' u^oo logU L « J - ' 

Substituting this estimate into (6.24) yields 

(6.32) limsup Ex, [Rk1{t^<k}] < -Ar? < 0. 

M— >oo iOg U 

Finally, substituting this last estimate and (6.12) into (6.5), we conclude that 

(6.33) lim inf ^E, [f:^ ] > -2^ + At?, 

u^oo logU 

as required. 

Next suppose that r = 1. Thus we are in the second scenario, where we 
would like to study the occupation measure Ou of {Vn} over an interval 
G{u) = [u^,u*] between the time T„ and K. By the Markov property, it is 
then sufficient to study a process {Vn} which begins at a level Vq > u and 
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terminates at time K. Suppose that Vq/u = v, for some constant v > 1. Then 
begin by writing (6.27) in a shghtly different form, namely, 

n oo oo 

(6.34) sup F„ - Von^^l - ^'' ^^^^^ ^' ■=^^i n ^J' 

1=1 1=1 j;=i+l 

provided that {Vq, . . . , Ki-i} is nonnegative. Now in the original measure, W 
converges to a proper random variable and hence for any 6 > 0,P (W > u ) — t- 
as u — )■ oo (c/. the proof of Lemma 5.2 in Collamore and Vidyashankar 
(2011)). Thus it is sufficient to study the occupation measure of the mul- 
tiplicative process VbHILi^* ^^ ^ region G'{u) = [n* ,u*], where [s',t'] is 
strictly contained in [s,t]. 

Let T„ denote the first time that the multiplicative random walk enters the 
set (— oo,u'' ]. Then by Heyde's (1966) a.s. version of the renewal theorem, 
we have for any e > that 



(6.35) P(T„>(^-^'-^)^°S" 



E[log^] 



1 as ti — )■ oo, 
u 



uniformly in w > 1. Similarly, if £„ denotes the last exit time of the multi- 
plicative random walk from the interval [u* ,oo), then by the strong law of 
large numbers, it follows that 



U as ti — 7- oo, 
u 



(636) P(z < (l-^^ + ^)l°g^ 

uniformly m v > 1. From the last two equations, we conclude that if ry > 
{t' - s' - 2e)/E[log^], then P (A^* > 77logn) -^ 1 as n -^ oo, where N* 
denotes the number of visits of the multiplicative random walk to the set 
G'{u). Consequently, 

(6.37) liminf— *— Eto(^)l >??, 

^ ' u^oo \ogu L M J - " 

yielding (6.33), as required. 

It remains to consider the case where {A) is violated. But mimicking the 
argument in Part I, we may set v^'^> = u + e/xj) and repeat the previous 
arguments to obtain that (6.33) holds, but with —2^ + Arj + eQi on the right- 
hand side, for some constant ©i (c/. (6.21)). Since e was arbitrary, the desired 
result follows. 

Finally, to complete the proof of theorem, note that if we do not have 
I' <^ ^j), as we have assumed throughout this proof, then by an application 
of the Radon-Nikodym theorem, 

h' = h'a + i^s, where Ua ^ /"D and u^ -L n^i- 



56 J. F. COLLAMORE ET AL. 

The proof can now be repeated, replacing everywhere v with Va- The mutuahy 
singular measure Vs plays no role, since -^ is zero on the support of this 
measure and thus, for example, the expectation on the right of (6.2) can be 
evaluated equally well with respect to Va as u. Moreover, with respect to 
the measure i>a, -^ = -^- The proof carries through in exactly the same 
manner as before, so we omit the details. □ 



7. Concluding remarks. In this paper, we developed and analyzed a 
dynamic importance sampling algorithm for rare event simulation of processes 
generated via SFPEs, encompassing a class of financial time series models 
and actuarial risk models. Extensions of these ideas to stratified importance 
sampling and general nonlinear recursions are currently under study. 
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